Analisis de datos omicos

1 Presentación

Este documento es una versión actualizada -en varios sentidos- de unos materiales de análisis de microparrays que escribimos allá por 2010, usando LateX y Sweave, junto con la profesora M. Carme Ruíz de Villa para la asignatura “Análisis de Microarrays” del posgrado de Bioinformática de la Universitat Oberta de Catalunya .

En estos años han cambiado muchas cosas. El posgrado es ahora un máster interuniversitario entre la UOC y la Universidad de Barcelona (UB). Muchos estadísticos i/o bioinformáticos -entre los que me incluyo- han cambiado de Latex/Sweave a RStudio/RMarkdown y Github. Mi compañera de escritura, Mamen, se ha jubilado, pero un nuevo compañero, Ricardo, se ha apuntado a la aventura de revivir estos materiales. Y aunque sigue sin estar claro si los microarrays también se jubilaran pronto, el campo de las ómicas ha crecido de forma explosiva. El “Análisis de datos ómicos”, título de este nuevo documento ya no se centra en los microarrays. También trata de RNA-seq, Proteómica, Metagenómica, Metabolómica o Epigenómica por citar sólo algunos de los más populares.

Es obvio que un curso de “Análisis de datos ómicos” no puede abarcar todas estas disciplinas, aunque también lo es que comparten muchos elementos comunes.

El plan de desarrollo para esta nueva edición, que podra seguirse en el repositorio de Github (), tendrá dos fases:

1. En una primera, actualizaremos la versión original a RMarkdown, eliminando aspectos desfasados y actualizando el código R a la versión actual de R/Bioconductor.
2. En una posterior miraré de añadir capítulos que introduzcan nuevas tecnologías y muestren como analizar los datos que estas generan. Si logro engañar a algunos colegas para que escriban un capítulo también los añadiré a los autores.

Previsiblemente este documento, con sus errores e imperfecciones se quedará como “open source” para uso de quien lo desee sin tener que pasar por el doloroso y anquilosante proceso de convertirlo en un libro, que, en este caso, empezaría a quedar obsoleto al minuto de estar acabado.

2 El proceso de análisis de microarrays

Este capítulo es una corta transición entre la primera parte del curso, en la que se han presentado los conceptos y herramientas básicos y la segunda en donde se presentan por separado y con mayor detalle los métodos de análisis de datos de microarrays.

Su objetivo por tanto es ofrecer una visión de conjunto que sirva de guía (“roadmap”) para los capítulos siguientes de forma que sin perder el detalle de cada uno de ellos tengamos conciencia de en que punto del proceso general nos encontramos.

El capítulo se estructura en dos partes. En la primera se presentan brevemente algunos de los problemas que típicamente se puede querer estudiar con microarrays u otras técnicas similares de análisis de datos de alto rendimiento. A continuación se presenta lo que se ha llamdo aquí el proceso de análisis de microarrays. Finalmente se intoducen algunos casos reales que, a modo de ejemplo se utilizarán en los capítulos siguientes.

2.1 Tipos de estudios

Los microarrays y otras tecnologías de alto rendimiento se han aplicado a multitud de investigaciones, de tipus muy diversos que van desde estudio del cancer (Alizadeh et al. (2000), Golub et al. (1999), van ’t Veer et al. (2002)) al de la germinación y la maduración del tomate (Moore et al. (2002)). A pesar de ello no resulta complicado clasificar los estudios realizados en algunos de los grandes bloques que se describen a continuación. La clasificación está basada en el excelente texto de Simon y colegas (Simon et al. (2003)) y aunque se origina en problemas de microarrays se puede aplicar fácilmente a estudios de genómica o ultrasecuenciación.

2.1.1 Comparación de grupos o “Class comparison”

El objetivo de los estudios comparativos es determinar si los perfiles de expresión génica difieren entre grupos previamente identificados. También se conoce estos estudios como de selección de genes diferencialmente expresados y son, sin duda los más habituales. Los grupos pueden representar una gran variedad de condiciones, desde distintos tejidos a distintos tratamientos o múltiples combinaciones de factores experimentales.

El análisis de este tipo de experimentos, que se describe en el capítulo sobre selección de genes diferencialmente expresados utiliza herramientas estadísticas como las pruebas de comparación de grupos paramétricas (t de Student) o no (test de Mann-Whitney) o diversos métodos de análisis de la varianza.

Entre los ejemplos de la sección @ref{c4examples} los casos @ref{celltypes}, @ref{estrogen} o @ref{CCL4} hacen referencia a estudios comparativos.

2.1.2 Predicción de clase o “Class prediction”

La predicción de clase puede confundirse con la selección de genes en tanto que disponemos de clases predefinidas pero su objetivo es distinto, ya que no pretende simplemente buscar genes cuya expresión sea distinta entre dichos grupos sino genes que puedan ser utilizados para identificar a que clase pertenece un “nuevo” individuo dado cuya clase es “a priori” desconocida. El proceso de predicción suele empezar con una selección de genes informativos, que pueden ser, o no, los mismos que se obtendrían si aplicáramos los métodos del apartado anterior, seguida de la construcción de un modelo de predicción y, lo que es más importante, de la verificación o validación de dicho modelo con unos datos nuevos independientes de los utilizados para el desarrollo del modelo.

Aunque el interés de la predicción de clase es muy alto se trata de un procedimiento mucho más complejo y con más posibilidades de error que la simple selección de genes diferencialmente expresados.

Entre los ejemplos de la sección @ref{c4examples} el caso @ref{golub} trata de un problema de predicción, a la vez que uno de descubrimiento de clases.

2.1.3 Descubrimiento de clases o “Class discovery”

Un problema distinto a los descritos se presenta cuando no se conoce las clases en que se agrupan los individuos. En este caso de lo que se trata es de encontrar grupos entre los datos que permitan reunir a los individuos más parecidos entre si y distintos de los de los demás grupos. Los métodos estadísticos que se emplearan en estos casos se conocen como análisis de clusters y no son tan complejos como los de predicción de clase aunque algunos aspectos como por ejemplo la definición del número de grupos no resulta tampoco sencillo.

Entre los ejemplos de la sección tanto el caso golub, en parte, como el @ref{breastTum} tratan problemas de descubrimiento de clases.

Una curiosidad del campo de la estadística es que el término clasificación aparece usado de forma indistinta para referirse a problemas de predicicón de clase o de descubrimiento de clase.

2.1.4 Otros tipos de estudios

Una vez identificados los principales tipos de estudios quedan muchos que no coinciden plenamente con ninguno de ellos. Sin entrar en detalles podemos señalar los estudios de evolución a lo largo del tiempo (“time course”), los de significación biológica (“Gene Enrichment Analysis”, “Gene Set Enrichment Analysis”, …) , los que buscan relaciones entre los genes (“network analysis” o “pathway analysis”). De momento con conocer e identificar los tres grandes bloques mencionados resultará más que suficiente.

2.2 Algunos ejemplos concretos

Una de las dificultades con que se encuentra la persona que comienza en el análisis de datos de microarrays es de donde obtener ejemplos concretos con los que prácticar las técnicas que está aprendiendo.

No es difícil encontrar datos de microarrays en internet por lo que se han seleccionado algunos conjuntos de datos interesantes para utilizarlos de ejemplo a lo largo del curso. Algunos de éstos son “populares” en el sentido de que han sido utilizados en diversas ocasiones y por lo tanto se encuentran bien documentados. Otros se han escogido simplemente porque ilustran bien algunas de las ideas que se desea exponer o por su accesibilidad.

Todos los datos corresponden a investigaciones publicadas por lo que no se describen exhaustivamente sinó que se expone brevemente el origen y objetivos del trabajo –incluyendo su clasificación según los grupos definidos en la sección anterior– y las características perinentes para el análisis como el tipo de microarrays, los grupos –si los hay– o como acceder a los datos.

2.2.1 Estudio de procesos regulados por citoquinas

Este conjunto de datos, que se denominará “celltypes”, corresponde a un estudio realizado por Chelvarajan y sus colegas (Chelvarajan et al. (2006)) que analizaron el efecto de la estimulación con lipopolisacáridos en la regulación por parte de citoquinas de ciertos procesos biológicos relacionados con la inflamación.

Este estudio es del tipo “class comparison” es decir su principal objetivo es la obtención de genes diferencialmente expresados entre dos o más condiciones.

Los datos se encuentran disponibles en la base de datos pública GEne Expression Omnibus mantenida por el National Institute of Health (NIH), pero pueden descargarse de la página de materiales del curso para garantizar su disponibilidad.

2.2.2 Clasificación molecular de la leucemia

A finales de los años 90, Todd Golub y sus colaboradores (Golub et al. (1999)) realizaron uno de los estudios más populares hasta el momento con datos de microarrays. En él utilizaron microarrays de oligonucleótidos para 6817 genes humanos para mirar de encontrar una forma de distinguir (clasificar) tumores de pacientes con leucemia linfoblástica aguda (ALL) de aquellos que sufrían de leucemia mieloide aguda (AML). Además se interesaba por la posibilidad de descubrir subgrupos de forma que pudieran definirse variantes de cada una de estas patologías a nivel molecular.

La diversidad de objetivos del estudio lleva a clasificarlo tanto entre los del tipo de predicción de clase como entre los que buscan descubrir nuievas clases o grupos en los datos.

Los datos de este estudio se encuentran disponibles en la web del instituto Broad, en donde se llevó a cabo (http://www.broadinstitute.org). También se encuentra disponible un paquete de R denominado ALL que permite utilizarlos directamente usando R y Bioconductor.

2.2.3 Efecto del estrogeno y el tiempo de administración {#estrogen]

Scholtens y colegas (Scholtens et al. (2004)) describen un estudio sobre el efecto de un tratamiento con estrógenos y del tiempo transcurrido desde el tratamiento en la expresión génica en pacientes de cáncer de mama. Los investigadores supusieron que los genes asociados con una respuesta temprana podrían considerarse dianas directas del tratamiento, mientras que los que tardaron más en hacerlo podrían considerarse objetivos secundarios correspondientes a dianas más alejadas en las vías metabólicas.

Estos datos han sido utilizados multitud de veces en los cursos de análisis de microarrays realizados por el proyecto Bioconductor y se encuentran disponibles en forma de paquete de R, el paquete estrogen. Una cracaterística importante de este paquetes es el hecho de que en vez de los datos procesados proporciona los datos “crudos”en forma de archivos .CEL de Affymetrix. Esto permite una mayor flexibilidad a la hora de reutilizarlos lo que explica su popularidad.

2.2.4 Efecto del CCL4 en la expresión génica

Holger y colegas de la empresa LGC Ltd. en Teddington, Inglaterra realizaron unos expeimentos con microarrays de dos colores en los que trataron hepatocitos de ratón con tetraclorido de carbono (CCL4) o con dimetilsulfóxido (DMSO). El tetraclorido de carbono fue ampliamente utilizado en productos de limpieza o refrigeración para el hogar hasta que se detectó que podía tener efectos tóxicos e incluso cancerígenos. El DMSO es un solvente similar, sin efectos tóxicos conocidos, que se utilizó como control negativo.

Los datos de este estudio no han sido publicados pero se encuentran disponibles en el paquete CCL4 de bioconductor. Su interés reside por un lado en que se trata de datos de microarrays de dos colores de la marca Agilent –en un momento en que la mayoría de estudios se realizan con datos de un color. Aparte de esto cabe resaltar el hecho de que el paquete incluye, de forma similar al anterior, los datos “crudos” en forma de archivos de tipo “Genepix” uno de los programas populares para escanear imágenes generadas con microarrays de dos colores.

Este estudio es también un estudio de comparación de clases cuyo objetivo principal es la selección de genes cuya expresión se asocia al tratamiento con CCL4 o DMSO.

2.2.5 Análisis de patrones en el ciclo celular

Los datos de este ejemplo denominado kidney son datos ya normalizados referidos a la expresión de los genes en distintos momentos del ciclo delular de la levadura e decir desde que concluye la división celular hasta que se inicia la siguiente.

Los datos puede descargarse desde la página del proyecto “Yeast Cell Cycle Project” (Proyecto de estudio del ciclo celular de la levadura) en la dirección:

http://genome-www.stanford.edu/cellcycle/data/rawdata/.

2.2.6 Recapitulación

La tabla resume la lista de ejemplos que se utilizan en este manual indicando el nombre con que nos referiremos de aquí en adelante a cada conjunto de datos as' como algunas de sus carcaterísticas.

La tabla siguiente resume los “datasets” utilizados en este manual. Aparte del nombre (arbitrario y “mnemotécnico”) se indica el tipo de microarrays, el número de muestras, y el tipo de problema para el que se utilizaron originalmente.

Nombre	Tipo (Modelo)	N. Muestras	Tipo de estudio
celltypes	Un color (Affy, Mouse 4302)	12	Comparativo
golub	Un color (Affy, HGU95A)	38	Clasificación
estrogen	Un color (Affy, HGU95A)	8	Comparativo
CCL4	Dos colores (Agilent, WG Rat Microarray)	16	Comparativo
breastTum	Un color (Affy, HGU95A)	49	Clasificación

2.3 El proceso de análisis de microarrays

Una vez descritos los tipos de análisis y algunos ejemplos podemos pasar a describir el proceso de análisis de microarrays que se resume brevemente en la figura @ref{c04analysisProcess}.

El análisis de microarrays, como la mayoría de análisis debe proceder de forma ordenada y siguiendo el método científico:

• La pregunta y su contexto nos servirán de guía para definir el Diseño experimental adecuado.
• El experimento se deberá realizar siguiendo las pautas decididas en el Diseño experimental y los datos obtenidos, que solemos denominar datos “crudos” o “raw data”, deberán someterse a los Controles de calidad adecuados antes de continuar con su análisis.
• Una vez decidido si la calidad de los datos es aceptable pasaremos a prepararlos para el análisis lo que puede incluir diversas formas de preprocesado, o transformaciones que a menudo se incluyen de forma general bajo el paraguas del término normalización, aunque, como veremos se trata de conceptos distintos.
• Los datos normalizados se utilizarán para los análisis estadísticos que hayamos decidido realizar durante el diseño experimental.
• Finalmente los resultados de los análisis serán la base para una interpretación biológica de los resultados del experimento.

Tal como ilustra la figura @ref(fig:c04analysisProcess>) el análisis de microarrays puede ser fácilmente visualizado como un proceso que empieza por una pregunta biológica y concluye con una interpretación de los resultados de los análisis que, de alguna forma, confiamos nos acerque un poco a la respuesta de la pregunta inicial.

Figure 2.1: Proceso de análisis de microarrays

El proceso descrito es básicamente una forma razonable de proceder en general. Los microarrays y otros datos genómicos son diferentes en su naturaleza de los datos clásicos alrededor de los que se han desarrollado la mayor parte de técnicas estadísticas. En consecuencia, en muchos casos ha sido necesario adaptar las técnicas existentes o desarrollar otras nuevas para adecuarse a las nuevas situaciones encontradas. Esto ha determinado que existan muchos métodos para cada una de los pasos descritos anteriormente lo que da lugar a una grandísima cantidad de posibilidades.

En la práctica lo que suele hacerse es optar por utilizar algunos de los métodos en los que hay un cierto consenso acerca de su calidad y utilidad para cada problema. Allison (Allison et al. (2006)) repasa los puntos principales de este consenso dando una lista de puntos a tener en cuenta en cualquier estudio que utilice microarrays. Imbeaud y Auffray (Imbeaud and Auffray (2005)) citan una lista de hasta 39 puntos que uno debe seguir en un experimento con microarrays para usar “buenas prácticas”.

Finalmente Zhu y otros (Zhu, Miecznikowski, and Halfon (2010)) utilizan un conjunto de arrays con valores de expresión conocidos para proponer los que, a su parecer, resultan los métodos más apropiados para cada etapa desde la corrección del backround hasta la selección de genes diferencialmente expresados. La figura 2.2 ilustra algunas de las opciones sugeridas por dichos autores.

Figure 2.2: Diseño de arrays

Los capítulos que siguen al presente proceden aproximadamente en el orden del proceso descrito en 2.1.

• Se empieza por tratar los principios del diseño de experimentos.
• A continuación se describen algunos métodos para el control de calidad, el preprocesado y la normalización de los datos.
• Se sigue con los métodos de selección de genes, que tratan aproximaciones básicas como el t-test y otras más sofisticadas como los modelos lineles para microarrays.
• La parte de selección de genes finaliza con una introducción a los métodos de análisis de la significación biológica de las listas de genes obtenidas de los procesos anteriores.
• Un último bloque aborda los métods de clasificación supervisada y no supervisada, que han tenido y siguen teniendo gran aplicabilidad en el análisis de todo tipo de datos ómicos.

3 Diseño de experimentos de microarrays

En este capítulo se examinan unas componentes clave del análisis de microarrays el diseño de experimentos que, no tan solo es crucial para una buena recogida de información, sino que marca todo el proceso desde el preprocesado al análisis final.

3.1 Fuentes de variabilidad

Los datos genómicos son muy variables. La figura 3.1, tomada de (???) ilustra algunas posibles fuentes de variabilidad, desde que se inicia el experimento hasta que se lee la información con el escáner. Sin tener que entrar en cómo influye cada una de estas fuents o cuan grande es su influencia lo que muestra esta figura es que este tipo de experimentos se enfrenta a múltiples fuentes de error, por lo que puede beneficiarse de un correcto diseño de experimentos.

Figure 3.1: Proceso de análisis de microarrays

3.2 Tipos de variabilidad

Habitualmente, en la mayor parte de situaciones experimentales, podemos distinguir entre variabilidad sistemática y variabilidad aleatoria.

La variabilidad sistemática es principalmente debida a procesos técnicos mientras que la variabilidad aleatoria es atribuible tanto a razones técnicas como biológicas. Podemos encontrar ejemplos de variabilidad sistemática en la extracción del ARN, marcaje o foto detección. La variabilidad aleatoria puede estar relacionada con muchos factores tales como la calidad del ADN o las características biológicas de las muestras.

La manera natural de manejar la variabilidad aleatoria es, por supuesto, la utilización de un diseño experimental apropiado y el apoyo para obtener conclusiones de unas herramientas estadísticas adecuadas. Los problemas relacionados con el c de los experimentos los discutiremos en esta sección y los relacionados con la aplicación de métodos estadísticos serán tratados en la sección métodos estadísticos.

Tradicionalmente, se estiman las correcciones de la variabilidad sistemática a partir de los datos, en lo que se llama genéricamente “calibrado”. En este contexto, hablaremos de “normalización”, que se tratará más adelante, en el capítulo dedicado al preprocesado de los datos.

3.3 Conceptos básicos de diseño de Experimentos

Empezaremos por definir conceptos que aparecen de forma usual cuando se plantea un diseño:

• Unidad experimental: Entidad física a la que se aplica un tratamiento, de forma independiente al resto de unidades. En cada unidad experimental se pueden realizar una medida o varias medidas, en este caso distinguiremos entre unidades experimentales y unidades observacionales.
• Factor: son las variables independientes que pueden influir en la variabilidad de la variable de interés.
  • Factor tratamiento: es un factor del que interesa conocer su influencia en la respuesta.
  • Factor bloque: es un factor en el que no se está interesado en conocer su influencia en la respuesta, pero se supone que esta existe y se quiere controlar para disminuir la variabilidad residual.
• Niveles : cada uno de los resultados de un factor. según sean elegidos por el experimentador o elegidos al azar de una amplia población se denominan factores de efectos fijos o factores de efectos aleatorios.
• Tratamiento : es una combinación especifica de los niveles de los factores en estudio. Son, por tanto, las condiciones experimentales que se desean comparar en el experimento. En un diseño con un único factor son los distintos niveles del factor y en un diseño con varios factores son las distintas combinaciones de niveles de los factores.
• Tamaño del Experimento: es el número total de observaciones recogidas en el diseño.

3.4 Principios de diseño de experimentos

Al planificar un experimento hay tres principios básicos que se deben tener siempre en cuenta: La replicación, el control local o “bloqueo” y la aleatorización.

Aplicados correctamente dichos principios garantizan diseños experimentales eficientes.

3.4.1 Replicación

La replicación o repetición de un experimento de forma idéntica en un número determinado de unidades, es la que permite la realización de un posterior análisis estadístico.

La utilización de replicas es importante para incrementar la precisión, obtener suficiente potencia en los tests y como base para los procedimientos de inferencia. Normalmente, distinguimos dos tipos de replicas en el análisis de microarrays:

• La replicación técnica se utiliza cuando estamos tratando réplicas del mismo material biológico. Puede ser tanto la replicación de spots en el mismo chip como la de diferentes alícuotas de la misma muestra hibridadas en diferentes microarrays.

• La replicación biológica se da cuando se toman medidas sobre muestras independientes, normalmente sobre individuos diferentes.

La replicación técnica permite la estimación del error a nivel de medida, mientras que la replicación biológica permite estimar la variabilidad a nivel de población.

Figure 3.2: Replicas biológicas vs réplicas técnicas

3.4.1.1 Potencia y tamaño de la muestra

Sorprendentemente, los primeros experimentos de microarrays utilizaban pocas o ninguna replica biológica. La principal explicación para este hecho, además de la falta de conocimiento estadístico, fue el alto coste de cada microarray.

En pocos años, la necesidad de las réplicas ha llegado a ser indiscutible, y al mismo tiempo, el coste de los chips ha disminuido considerablemente. Actualmente, es normal la utilización de, al menos, tres a cinco replicas por condición experimental, aunque a este consenso se ha llegado más por razones empíricas que por la disponibilidad de modelos apropiados para el análisis de la potencia y tamaño de la muestra.

Recientemente, se ha producido una importante afluencia de artículos describiendo métodos para el análisis de la potencia y tamaño de la muestra. A pesar de su variedad, ningún método aparece como candidato claro para ser utilizado en situaciones prácticas. Esto se debe, probablemente, a la complejidad de los datos de microarrays, básicamente por que los genes no son independientes. Por tanto, estas estructuras de correlación existen en los datos, pero la mayor parte de las dependencias son desconocidas por lo que la estimación de estas estructuras es muy complicada.

Como indicaba Allison Allison et al. (2006), aunque no hay consenso sobre que procedimiento es mejor para determinar el tamaño de las muestras, sí que lo hay sobre la conveniencia de realizar el análisis de la potencia, y, por supuesto, sobre el hecho de que un mayor número de replicas generalmente proporcionan una mayor potencia.

3.4.1.2 Pooling

En el contexto de los microarrays, llamamos “pooling” a la combinación del RNA de diferentes casos en una única muestra.

Hay dos razones a favor de ello, una, es que, a veces, no hay suficiente ARN disponible y esta es la única forma de conseguir suficiente material para construir los arrays, otra, más controvertida, es la creencia que la variabilidad entre arrays puede reducirse por “pooling”. La justificación es que combinar muestras equivale a promediar expresiones, y como ya se conoce, el promedio es menos variable que los valores individuales. A pesar de la debilidad de este argumento, es verdad que en ciertas situaciones el pooling puede ser apropiado y, recientemente, muchos estadísticos han dedicado sus esfuerzos para tratar de responder la pregunta “to pool or not to pool” (Kerr and Churchill 2001). Por ejemplo, si la variabilidad biológica está altamente relacionada con el error en las medidas, y las muestras biológicas tienen un coste mínimo en comparación con el de los arrays, una apropiada estrategia de “pooling” puede ser claramente eficiente en costes.

De todos modos, el “pooling” no se debería usar en cualquier tipo de estudios. Si el objetivo es comparar expresiones medias (ver más adelante " comparación de clases"), puede funcionar adecuadamente, pero se debería claramente evitar cuando el objetivo del experimento es construir predictores que se basen en características individuales.

3.4.2 Aleatorización

Se entiende por aleatorizar la asignación de todos los factores al azar a las unidades experimentales. Co ello se consigue disminuir el efecto de los factores no controlados por el experimentador en el diseño experimental y que podrían influir en los resultados.

Las ventajas de aleatorizar los factores no controlados son:

• Transforma la variabilidad sistemática no planificada en variabilidad no planificada o ruido aleatorio. Dicho de otra forma, aleatorizar previene contra la introducción de sesgos en el experimento.
• Evita la dependencia entre observaciones al aleatorizar los instantes de recogida muestral.
• Valida muchos de los procedimientos estadísticos más comunes.

3.4.3 Bloqueo o control local

Hace referencia a dividir o particionar las unidades experimentales en grupos llamados bloques de modo que las observaciones realizadas en cada bloque se realicen bajo condiciones experimentales lo más parecidas posibles. A diferencia de lo que ocurre con los factores tratamiento, el experimentador no está interesado en investigar las posibles diferencias de la respuesta entre los niveles de los factores bloque.

Bloquear es una buena estrategia siempre y cuando sea posible dividir las unidades experimentales en grupos de unidades similares. La ventaja de bloquear un factor que se supone que tiene una clara influencia en la respuesta, pero en el que no se está interesado, es que convierte la variabilidad sistemática no planificada en variabilidad sistemática planificada.

3.5 Diseños experimentales para microarrays de dos colores

En arrays de dos colores, se aplican dos condiciones experimentales a cada array. Esto permite la estimación del efecto del array, como el efecto de bloqueo. En Affymetrix u otro array de un canal, cada condición se aplica a un chip separado, imposibilitando la estimación del efecto de los microarrays, lo cual, por otra parte, se considera que tiene una relación muy pequeña en el tratamiento de los efectos debido al proceso industrial utilizado para fabricar estos chips.

Como consecuencia de lo anterior, los experimentos que usan arrays de un canal se consideran “estándar”, por lo que se les pueden aplicar los conceptos y técnicas tradicionales de diseño de experimentos .

Los arrays de dos canales presentan una situación más complicada. Por una parte, los “dos colores” no son simétricos, es decir, con la misma cantidad de material, un array hibridado con uno u otro color sea Cy5 o Cy3, emitirá señales con diferente intensidad. La forma normal de manejar este problema es dye-swapping que consiste en utilizar para una misma comparación dos arrays con dyes cambiados, es decir, si en el primer array la muestra 1 se marca con Cy3 y la muestra 2 con Cy5, con las muestras del segundo array se hace al revés (ver figura 3.3.

Figure 3.3: Representación simplificada de dos diseños

Por otra parte, el hecho de que solo se puedan aplicar dos condiciones a cada array complica el diseño, ya sea porque normalmente hay más de dos condiciones, o porque no es recomendable hibridar directamente dos muestras en un array, creando pares artificiales.

El problema de la asignación eficiente de muestras a microarrays, dado un numero de condiciones a comparar y un número fijo de arrays disponibles, ha sido estudiado de forma exhaustiva.

El diseño utilizado más comúnmente en la comunidad biológica es el diseño de referencia en el que cada condición de interés se compara con muestras tomadas de alguna referencia estándar común a todos las arrays, (ver figura 3.4, imagen (a)).

Figure 3.4: (a) diseño de referencia. (b) diseño en loop.

Los diseños de referencia permiten hacer comparaciones indirectas entre condiciones de interés. La crítica más importante a esta aproximación es que el 50 de las fuentes de hibridación se utilizan para producir la señal del grupo control, de un interés no intrínseco para los biólogos. En contraste, un diseño en loop compara dos condiciones a través de otra cadena de condiciones, por lo que elimina la necesidad de una muestra de referencia.

4 Exploración de los datos, control de calidad y preprocesado

4.1 Introducción

Los estudios realizados con microarrays, sea cual sea la tecnología en que se basan, tienen una característica común: generan grandes cantidades de datos a través de una serie de procesos, que hacen que su significado no siempre sea completamente intuitivo.

Como en todo tipo de análisis, antes de empezar a trabajar con los datos, debemos de asegurarnos de que éstos son fiables y completos y de que se encuentran en la escala apropiada para proporcionar la información que pretendemos obtener de ellos.

En el caso de los microarrays solemos distinguir dos fases previas al análisis de los datos:

1. Exploración y control de calidad Mediante gráficos y/o resúmenes numéricos estudiamos la estructura de los datos con el fin de decidir si parecen correctos o presentan anomalías que deben ser corregidas.

2. Preprocesado Incluso si son correctos, los datos “crudos” no sirven para el análisis sino que necesitan preprocesados, lo que puede interpretarse como uno o más de los procesos siguientes:

• “Limpiados”, para eliminar la parte de la señal no atribuible a la expresión, llamada “background” o ruído.
• “Normalizados” para hacerlos comparables entre muestras y eliminar posibles sesgos técnicos.
• Resumidos o “sumarizados” de forma que se tenga un sólo valor por gen.
• “Transformados” de forma que la escala sea razonable y facilite el análisis.

A menudo -por un cierto abuso de lenguaje- se denomina “normalización” al preprocesado descrito en las etapas anteriores.

4.1.1 Nivel de análisis y tipo de microarray

Desde la generalización del uso de los microarrays se han desarrollado muchas formas para visualizar los datos y decidir acerca de su calidad. Algunas trabajan sobre los datos obtenidos del escáner, otras lo hacen con los datos normalizados. Algunas sirven tanto para arrays de dos colores como arrays de un color. Otras son específicas de la tecnología. Con el fin de organizar esta multitud de opciones podemos diferenciar:

• Datos de bajo nivel Son los datos proporcionados por el escáner contenidos por ejemplo en archivos .gpr (dos colores) o .cel (un color). Estos últimos son binarios por lo que no es posible ni tan sólo visualizarlos sin programas específicos.
• Datos de alto nivel Son los datos resultantes del (pre)procesado de los datos de bajo nivel. Básicamente se corresponden con los datos de expresión ya resumidos (“sumarizados”), normalizados o no.

La exploración y el control de calidad pueden basarse en datos de bajo o de alto nivel. En cada caso se pueden aplicar ya sean tácnicas generales de visualización de datos, como histogramas, diagramas de caja o visualización en dimensiones reducidas (PCA u otras), o bien técnicas “ad-hoc” para cada tipo de datos como el MA-Plot u otras que se discuten a continuación.

4.1.2 Datos de partida

La estructura de los datos de microarrays de un color o de dos colores difiere considerablemente, tanto a nivel físico (los chips y las imágenes que de ellos se obtiene) como informático, es decir en la forma en que se representan.

4.1.2.1 Arrays de dos colores (o de “cDNA”)

Tradicionalmente los arrays de dos colores o de cDNA se realizaban de forma menos automatizada que los de un color o de Affymetrix. Esto implica que, tras obtener la imagen, el escaneado del archivo “.TIFF” resultante pudiera ser llevado a cabo mediante un software independiente como Genepix, un programa que convierte las imágenes en números y genera un archivo de información (con extensión “.gpr”) a partir del cual pueden calcularse las expresiones relativas, así como valores de calidad para cada “spot” o punto escaneado en la imagen.

Para cada imagen (o sea para cada microarray) hay un archivo “.gpr” que contiene una fila por gen y varias columnas con distintos valores, por ejemplo:

• la intensidad para cada canal,
• valores resumen de las intensidades y
• resultados de controles de calidad (“FLAGS”).

La tabla siguiente muestra lo que serían las primeras filas y columnas de un archivo “.gpr” obtenido mediante el programa Genepix.

Gen	Señal	Fondo	Señal	Fondo	“Flag”	Otros
	Cy5 (R)	Cy5 (Rb)	Cy3 (G)	Cy3 (Gb)
gen-1	3.7547	1.8128	5.0672	1.8496	1	…
gen-2	0.8331	0.9175	1.1536	0.9995	0	…
gen-3	9.8254	0.2781	0.6921	0.5430	1	…
gen-4	9.1539	0.1918	3.8290	0.0014	0	…
gen-5	4.8603	0.2377	0.5338	0.3335	0	…
gen-6	7.8567	1.3941	1.3050	1.4876	1	…
gen-7	2.7619	0.8108	8.4916	1.6518	0	…
gen-8	3.6618	0.1918	0.3770	1.1842	0	…
gen-9	4.4300	0.5888	2.8161	0.8707	1	…
…	…	…	…	…	…	…

Los valores de intensidad se convierten en una única matriz de expresión que contiene una columna por chip con los valores de intensidad relativa obtenidas por ejemplo con una “sumarización” del tipo: \[ \log\frac{R-Rb}{G-Gb}, \] y una fila por gen (mismas filas que archivos “.gpr”).

La tabla siguiente muestra lo que sería una matriz de expresión derivada de cuatro archivos “.gpr” como los de la tabla anterior.

gen	Array-1	Array-2	Array-3	Array-4
gen-1	1.65695	-0.10820	1.69515	8.25137
gen-2	-1.82305	0.11350	1.58807	0.95676
gen-3	0.01561	0.47682	-27.88036	-1.39078
gen-4	0.42709	4.29319	-4.31366	30.96866
gen-5	0.04332	0.24126	-10.27675	0.26680
gen-6	-0.02825	0.68408	2.04163	0.99554
gen-7	3.50549	-8.04635	0.18286	0.22647
gen-8	-0.23261	-1.08477	0.19582	1.11561
gen-9	0.50643	1.52147	0.99092	0.23441
…	…	…	…	…

4.1.2.2 Arrays de un color (Affymetrix)

El resultado de escanear la imagen de un array de affymetrix es un archivo de extensión “.CEL” que, a diferencia de los arrays de dos colores, está en formato binario es decir que solo puede ser leído con programas específicos para ello.

De forma similar a los arrays de dos colores, existe un archivo “.CEL” por cada microarray.

A partir de las intensidades de los archivos “.CEL” se genera la matriz de expresión, que contiene una columna por chip con los valores de intensidad absoluta, y una fila por grupo de sondas. En el caso de arrays de affymetrix existe una gran variedad de algoritmos de “sumarización” y, según cual se utilice, se obtendrá unos u otros valores de expresión. Ahora bien, éstos serán siempre medidas absolutas, es decir independientes del resto de muestras y en una escala arbitraria.

La tabla siguiente muestra las primeras filas de una matriz de expresión sumarizada correspondiente a los primeros genes de uno de los casos resueltos.

ID_REF	GSM188013	GSM188014	GSM188016	GSM188018
1007_s_at	15630.2	17048.8	13667.5	15138.8
1053_at	3614.4	3563.22	2604.65	1945.71
117_at	1032.67	1164.15	510.692	5061.2
121_at	5917.8	6826.67	4562.44	5870.13
1255_g_at	224.525	395.025	207.087	164.835
…	…	…	…	…

4.1.2.3 Datos de ejemplo

Los ejemplos de este capítulo se basarán en dos de los conjuntos de datos descritos en el capítulo 2.

• Por un lado, trabajaremos con los datos del conjunto referidos al efecto del tratamiento con estrógenos en cancer de mama descrito en @ref{estrogen}. La exploración se basará en los datos normalizados pero también en los datos “crudos”, de bajo nivel, contenidos en los archivos “.cel”. Estos datos pueden obtenerse directamente del paquete de Bioconductor . Un aspecto interesante de este conjunto de datos es que, junto con 8 muestras “buenas” se proporciona un array defectuoso, denominado “” que facilita el ver como aparecen ciertos gráficos cuando hay problemas.
• Por otro lado, usaremos los datos del conjunto , relativos al efecto del tratamiento con CCL4 en la expresión de los hepatocitos. Los datos resultan accesibles al instalar el paquete .

if(!require(BiocManager)) install.packages("BiocManager")
if (!(require(CCl4))){
 BiocManager::install("CCl4")
}
if (!(require(estrogen))){
 BiocManager::install("estrogen")
}
if (!(require(affy))){
 BiocManager::install("affy")
}
if (!(require(affyPLM))){
 BiocManager::install("affyPLM")
}

Aunque, en la actualidad (año 2021) el paquete más utilizado para arrays de este tipo es el paquete oligo este ejemplo utiliza el paquete affy.

library(estrogen)
library(affy)
affyPath <- system.file("extdata", package = "estrogen")
adfAffy = read.AnnotatedDataFrame("phenoData.txt", sep="",  path=affyPath)
affyTargets = pData(adfAffy)
affyTargets$filename = file.path(affyPath, row.names(affyTargets))
affyRaw <- read.affybatch(affyTargets$filename, phenoData=adfAffy)
# show(affyRaw)
actualPath <- getwd()
setwd(affyPath)
allAffyRaw <- ReadAffy()
setwd(actualPath)

El resultado de la lectura es un objeto “affyraw” de clase “affyBatch”:

class(affyRaw)

## [1] "AffyBatch"
## attr(,"package")
## [1] "affy"

print(affyRaw)

## AffyBatch object
## size of arrays=640x640 features (22 kb)
## cdf=HG_U95Av2 (12625 affyids)
## number of samples=8
## number of genes=12625
## annotation=hgu95av2
## notes=

El paquete limma es conocido como paquete para la selección de genes diferencialmente expresados pero también incluye algunas funciones para la lectura y preprocesado de arrays de dos colores.

library("limma")
library("CCl4")
dataPath = system.file("extdata", package="CCl4")
adf = read.AnnotatedDataFrame("samplesInfo.txt", 
    path=dataPath)
#adf
targets = pData(adf)
targets$FileName = row.names(targets)
RG <- read.maimages(targets, path=dataPath, source="genepix")
attach(RG$targets)
newNames <-paste(substr(Cy3,1,3),substr(Cy5,1,3),substr(FileName,10,12), sep="")
colnames(RG$R)<-colnames(RG$G)<-colnames(RG$Rb)<-colnames(RG$Gb)<-rownames(RG$targets)<- newNames
# show(RG)

El resultado de la lectura es un objeto de classe “RGlist”.

class(RG)

## [1] "RGList"
## attr(,"package")
## [1] "limma"

print(RG)

## An object of class "RGList"
## $R
##      DMSCCl319 DMSCCl320 DMSCCl321 DMSCCl329 CClDMS330 CClDMS331 CClDMS332
## [1,]        47        46        40        59        56        40        45
## [2,]      1095        46        39        54        56        40        47
## [3,]        47        46        41        55        52        39        46
## [4,]        46        45        41        56        54        43        47
## [5,]        55        53        45        61        71        47        63
##      DMSCCl333 CClDMS379 CClDMS380 CClDMS381 DMSCCl382 DMSCCl384 DMSCCl389
## [1,]        70        46        65        43        59        46        63
## [2,]        69        45        61        43        62        48        68
## [3,]        81        44        57        41        58        50        63
## [4,]        78        44        60        42        59        48        62
## [5,]        91        51        83        56        71        52        69
##      DMSCCl390 CClDMS391 CClDMS393 CClDMS394
## [1,]        77        40        91        32
## [2,]        71        41        85        32
## [3,]        70        40        39        32
## [4,]        81        41        51        32
## [5,]        81        57        95        37
## 44285 more rows ...
## 
## $G
##      DMSCCl319 DMSCCl320 DMSCCl321 DMSCCl329 CClDMS330 CClDMS331 CClDMS332
## [1,]        77        78        45        45        87        65       128
## [2,]      4558        79        44        47        81        66       131
## [3,]        79        75        46        47        84        66       122
## [4,]        78        82        44        71        86        63       135
## [5,]        93        95        92        57        86        66       132
##      DMSCCl333 CClDMS379 CClDMS380 CClDMS381 DMSCCl382 DMSCCl384 DMSCCl389
## [1,]       107        57        90        99        84        56        68
## [2,]       110        54        89        96       168        54        72
## [3,]       114        54        84        93        78        55        71
## [4,]       113       118        84        97        83        56        75
## [5,]       133        59        95        99       102        64        84
##      DMSCCl390 CClDMS391 CClDMS393 CClDMS394
## [1,]        71        74        83        59
## [2,]        71        75        81        60
## [3,]        73        78        62        59
## [4,]        79        77        70        61
## [5,]        99        83        92        63
## 44285 more rows ...
## 
## $Rb
##      DMSCCl319 DMSCCl320 DMSCCl321 DMSCCl329 CClDMS330 CClDMS331 CClDMS332
## [1,]        41        39        35        47        44        36        40
## [2,]        43        39        35        46        44        36        39
## [3,]        41        40        35        48        44        35        39
## [4,]        40        40        35        49        43        35        39
## [5,]        41        40        36        48        44        35        39
##      DMSCCl333 CClDMS379 CClDMS380 CClDMS381 DMSCCl382 DMSCCl384 DMSCCl389
## [1,]        55        39        49        37        46        42        50
## [2,]        55        37        50        37        46        41        52
## [3,]        59        39        47        37        46        41        52
## [4,]        63        39        48        36        46        41        52
## [5,]        60        39        47        37        46        41        52
##      DMSCCl390 CClDMS391 CClDMS393 CClDMS394
## [1,]        60        36        34        30
## [2,]        59        36        35        31
## [3,]        55        35        35        30
## [4,]        58        36        34        30
## [5,]        56        36        35        30
## 44285 more rows ...
## 
## $Gb
##      DMSCCl319 DMSCCl320 DMSCCl321 DMSCCl329 CClDMS330 CClDMS331 CClDMS332
## [1,]        44        44        32        34        40        39        48
## [2,]        46        43        32        35        41        38        48
## [3,]        43        43        32        35        40        38        49
## [4,]        43        43        32        35        41        38        49
## [5,]        43        42        32        35        41        38        48
##      DMSCCl333 CClDMS379 CClDMS380 CClDMS381 DMSCCl382 DMSCCl384 DMSCCl389
## [1,]        47        35        50        50        45        37        38
## [2,]        48        35        50        49        45        37        38
## [3,]        48        35        51        49        45        37        38
## [4,]        49        35        51        49        46        37        37
## [5,]        49        35        50        49        45        37        38
##      DMSCCl390 CClDMS391 CClDMS393 CClDMS394
## [1,]        38        40        37        36
## [2,]        38        40        37        36
## [3,]        39        40        37        36
## [4,]        38        40        37        36
## [5,]        39        40        37        36
## 44285 more rows ...
## 
## $targets
##            Cy3  Cy5 RIN.Cy3 RIN.Cy5                      FileName
## DMSCCl319 DMSO CCl4     9.0     9.7 251316214319_auto_479-628.gpr
## DMSCCl320 DMSO CCl4     9.0     5.0 251316214320_auto_478-629.gpr
## DMSCCl321 DMSO CCl4     9.0     2.5 251316214321_auto_410-592.gpr
## DMSCCl329 DMSO CCl4     9.0     2.5 251316214329_auto_429-673.gpr
## CClDMS330 CCl4 DMSO     9.7     9.0 251316214330_auto_457-658.gpr
## 13 more rows ...
## 
## $genes
##   Block Row Column           ID            Name
## 1     1   1      1 BrightCorner    BrightCorner
## 2     1   1      2 BrightCorner    BrightCorner
## 3     1   1      3    (-)3xSLv1 NegativeControl
## 4     1   1      4 A_44_P317301        AW523361
## 5     1   1      5 A_44_P386163    NM_001007719
## 44285 more rows ...
## 
## $source
## [1] "genepix"
## 
## $printer
## $ngrid.r
## [1] 1
## 
## $ngrid.c
## [1] 1
## 
## $nspot.r
## [1] 430
## 
## $nspot.c
## [1] 103

4.2 Exploración y control de calidad

Los gráficos son útiles para comprobar la calidad de los datos de microarrays, obtener información sobre cómo se deben preprocesar los datos y comprobar, finalmente, que el preprocesado se haya realizado correctamente.

Siguiendo el esquema presentado en la tabla @ref{tab:tablaDiagnosticos} se presentan a continuación los distintos gráficos utilizados con una breve descripción de lo que representa cada uno y como interpretarlos adecuadamente. El código para generarlos se presenta al final del capítulo como un apéndice.

4.2.1 Control de calidad con gráficos estadísticos generales

4.2.1.1 Histogramas y gráficos de densidad

Estos gráficos permite hacerse una idea de si las distribuciones de los distintos arrays son similares en forma y posición. La figura ?? muestra los histogramas correspondientes a los 9 arrays del conjunto de datos .

affySampleNames <- rownames(pData(allAffyRaw))
affyColores <- c(1,2,2,3,3,4,4,8,8)
affyLineas <- c(1,2,2,2,2,3,3,3,3)
hist(allAffyRaw, main="Signal distribution", col=affyColores, lty=affyLineas)
legend (x="topright", legend=affySampleNames , col=affyColores, lty=affyLineas, cex=0.7)

4.2.1.2 Diagramas de caja o “boxplots”

Como los histogramas los diagramas de caja –basados en los distintos cuantiles de las valores– dan una idea de la distribución de las intensidades. La figura ?? muestra los diagramas de caja correspondientes a los 9 arrays del conjunto de datos .

boxplot(allAffyRaw, main="Signal distribution", col=affyColores, las=2)

4.2.1.3 Gráficos de componentes principales

El análisis de componentes principales puede servir para detectar si las muestras se agrupan de forma “natural” es decir con otras muestras provenientes del mismo grupo o si no hay correspondencia clara entre grupos experimentales y proximidad en este gráfico. Cuando esto sucede no significa necesariamente que haya un problema pero puede ser indicativo de efectos técnicos -como el conocido efecto “batch”- que podría ser necesario corregir.

plotPCA <- function ( X, labels=NULL, colors=NULL, dataDesc="", scale=FALSE)
{
  pcX<-prcomp(t(X), scale=scale) # o prcomp(t(X))
  loads<- round(pcX$sdev^2/sum(pcX$sdev^2)*100,1)
  xlab<-c(paste("PC1",loads[1],"%"))
  ylab<-c(paste("PC2",loads[2],"%"))
  if (is.null(colors)) colors=1
  plot(pcX$x[,1:2],xlab=xlab,ylab=ylab, col=colors, 
       xlim=c(min(pcX$x[,1])-10, max(pcX$x[,1])+10),
       ylim=c(min(pcX$x[,2])-10, max(pcX$x[,2])+10),
       )
  text(pcX$x[,1],pcX$x[,2], labels, pos=3, cex=0.8)
  title(paste("Plot of first 2 PCs for expressions in", dataDesc, sep=" "), cex=0.8)
}

if (!file.exists("allAffyNorm")){
  allAffyNorm<- rma(allAffyRaw)
  affyNorm <- rma(affyRaw)
  save(allAffyNorm, affyNorm, file="affyNorm.Rda")
}else{
  load(file="affyNorm.Rda")
}

La figura ?? muestra dos diagramas de compnentes principales realizados a partir de los datos normalizados del conjunto de datos . El gráfico de la parte superior que incluye el array defectuoso ilustra que la principal fuente de variabilidad es la diferencia de este array con el resto. Cuando se repite el análisis omitiendo esta muestra puede verse como la principal fuente de variación (eje X) se asocia con el tiempo de exposición (alto a la derecha, bajo (10h) a la izquierda, mientras que la segunda fuente de variación se asocia con la exposición a los estrógenos (alto arriba, bajo abajo).

El gráfico de la parte inferior en el que este array se ha eliminado antes de calcular las componentes principales muestra dos agrupaciones claras de arrays “derecha e izquierda” y “arriba y abajo” de la figura. El hecho de que los grupos no se correspondan con los tratamientos sugiere que puede haber alguna fuente adicional de variación que ha de tenerse en cuenta.

opt <- par(mfrow=c(2,1))
plotPCA(exprs(allAffyNorm), labels=affySampleNames, dataDesc="PCA for all arrays\nincludes defective sample")
plotPCA(exprs(affyNorm), labels=colnames(exprs(affyNorm)), dataDesc="PCA for all arrays")

par(opt)

4.2.1.4 Imagen del chip

Otra forma de ver si las muestras se agrupan según los grupos experimentales, o mediante otros criterios es usando un cluster jerárquico que realiza una agrupación básica de las muestras por grado de similaridad según la distancia que se utilice.

Como en el caso de las componentes principales si las muestras se agrupan según las condiciones experimentales es una buena señal pero si no es así puede deberse a la presencia de otra fuente de variación o bien al hecho de que se trata de un gráfico basado en todo los datos y las condiciones experimentales pueden haber afectado un pequeño número de genes.

La figura 4.1 muestra como se agrupan los datos del conjunto en base a un cluster jerárquico. Como en el caso de las componentes principales tras eliminar el array defectuoso las muestras se separan, primero por el tiempo de exposición y luego por niveles de estrógeno suministrado.

clust.euclid.average <- hclust(dist(t(exprs(affyNorm))),method="average")
plot(clust.euclid.average, labels=colnames(exprs(affyNorm)), main="Hierarchical clustering of samples",  hang=-1, cex=0.7)

Figure 4.1: Un cluster jerárquico sirve para determinar si las muestras se agrupan de forma natural según los grupos experimentales o si lo hacen de otra forma

4.2.2 Gráficos de diagnóstico para microarrays de dos colores

El diagnóstico de arrays de dos canales se basa principalmente en la imagen y en diferentes tipos de gráficos.

4.2.2.1 Diagramas de dispersión y “MA-plots”

La normalización discutida en este mismo capítulo es un punto clave en el proceso de análisis de microarrays y se ha dedicado un gran esfuerzo a desarrollar y probar diferentes métodos (Quackenbush (2002),Yang et al. (2002)). Una razón para ello es que hay diferentes artefactos técnicos que deben ser corregidos para poder ser utilizados, y no cualquier método puede funcionar con todos ellos.

En general, los métodos de normalización se basan en el siguiente principio: .

Esto da una idea de como debería ser un gráfico de intensidades. Por ejemplo, si no hubiese artefactos técnicos, en arrays de dos canales, una gráfica de dispersión de intensidad del rojo frente al verde debería dejar la mayor parte de los puntos alrededor de una diagonal. Cualquier desviación de esta situación debería ser atribuible a razones técnicas, no biológicas, y por tanto, debería ser eliminada. Esto ha conducido a un método de normalización muy popular consistente en estimar la transformación a aplicar, como una función de las intensidades utilizando el método en la representación transformada de la gráfica de dispersión conocida como el .

La figura 4.2 muestra, en su parte izquierda (a) un gráfico de dispersión del canal rojo frente al verde en un array de dos colores. El hecho de que los datos no estén centrados alrededor de la diagonal sugiere la necesidad de normalización.

Una representación muy popular que ayuda a visualizar mejor esta asimetría es lo que se conoce como “”, que aparece en la parte derecha (b) de la figura4.2. Geométricamente representa una rotación del gráfico de dispersión, en la que el significado de los nuevos ejes es:

• \(A=\displaystyle \frac{1}{2}(\log_2 (R*G))\): El logaritmo de la intensidad media de los dos canales,
• \(M=\log_2 \displaystyle \frac{R}{G}\): El logaritmo de la expresión relativa entre ambos canales (normalmente conocido como “log–ratio”).

Figure 4.2: (a) Gráfico de un canal frente al otro (b) MA-plot (intensidad frente log-ratio)

4.2.2.2 Imagen del array

La imagen del chip (véase la figura ??, izquierda) ofrece una visión rápida de la calidad del array, proporcionando información acerca del balance del color, la uniformidad en la hibridación y en los , de si el background es mayor del normal y dela existencia de artefactos como el polvo o pequeñas marcas (rasguños).

4.2.2.3 Histogramas de señales y de la relación señal–ruído

Estos gráficos (véase la figura 4.3, derecha) son útiles para detectar posibles anormalidades o un background excesivamente alto .

Figure 4.3: La relación señal ruido sirve para detectar posibles anormalidades o un background excesivamente alto como medida de calidad

4.2.2.4 Boxplots

Un gráfico muy utilizado es el diagrama de cajas o “boxplot” múltiple con una caja por cada chip. Del alineamiento (o falta de él) y la semejanza (o disparidad) entre las cajas, se deduce si hace falta, o no, normalizar entre arrays.

En el caso de arrays de dos colores pueden utilizarse diagrams de cajas “dentro de arrays” (entre distintos sectores del mismo chip) y “entre arrays”.

4.2.3 Gráficos de diagnóstico para microarrays de un color

4.2.3.1 Imagen del chip

Los arrays de affymetrix contienen millones de sondas por lo que no pueden examinarse a simple vista. A pesar de ello hay diversas formas de obtener una imagen que, en caso de presentar irregularidades pueden indicar algún tipo de problemas como burbujas, arañazos, etc. La figura @ref(c06plotAffy} muestra algunas de las imágenes

La imagen del array de Affymetrix sólo es útil para evidenciar grandes problemas como burbujas, arañazos, etc. En este ejemplo los dos arrays de la izquierda se considerarían aceptables y los de la derecha defectuosos.

Figure 4.4: Imágenes de cuatro microarrays de Affymetrix

4.2.3.2 Gráfico “M-A”

En los chips de dos colores el MA–plot se utliza para comparar los dos canales en cada array (rojo y verde). En cambio, en los chips de Affymetrics, en que sólo hay un canal en cada array, la única forma de definir M (el log ratio) es a partir de la comparación entre pares de de valores, ya sea los arrays dos a dos o bien cada array respecto un valor de referencia que puede ser la mediana, punto a punto, de todos los arrays (véase por ejemplo la figure ).

Figure 4.5: En los chips de Affymetrix la única forma de definir M (el log ratio) es comparar entre diferentes arrays

\(M=\log_2(I_1) - \log_2(I_2)\): log ratio \(A=\displaystyle \frac{1}{2}(\log_2 (I_1)+\log_2(I_2))\): log de intensidades Donde \(I_1\) es la intensidad del array de estudio, e \(I_2\) es la intensidad media de arrays. Por lo general, se espera que la distribución en el gráfico se concentre a lo largo del eje M = 0.

4.2.3.3 Modelos de bajo nivel (“Probe-Level-Models” o PLM)

Los modelos de bajo nivel (“Probe-Level-Models” o PLM) ajustan a los valores de intensidad –a nivel de sondas, no de valores totalizados de gen– un modelo explicativo. Los valores estimados por este modelo se comparan con los valores reales y se obtienen los errores o “residuos” del ajuste. El análisis de dichos residuos procede de forma similar a lo que se realiza al analizar un modelo de regresión: Si los errores no presentan ningún patrón especial supondremos que el modelo se ajusta relativamente bien. Si, en cambio, observamos desviaciones de esta presunta aleatoriedad querrá decir que el modelo no explica bien las observaciones, lo cual se atribuirá a la existencia de algún problema con los datos.

Con los valores ajustados del modelo se calculan dos medidas:

• La expresión relativa en escala logarítmica " Relative Log Expression" (RLE) es una medida estandarizada de la expresión. No es de gran utilidad pero debería presentar una distribución similar en todos los arrays.
• El error no estandarizado y normalizado o “NUSE” es el más informativo ya que representa la distribucián de los residuos a la que hacíamos referencia más arriba. Si un array es problemático la caja correspondiente en el boxplot aparece desplazada hacia arriba o abajo de las demás.

stopifnot(require(affyPLM))
Pset<- fitPLM(allAffyRaw)

La figura ?? muestra los gráficos RLE y NUSE para el conjunto de datos estrogen. En ambos gráficos puede verse como el array defectuoso “” queda claramente diferenciado del resto.

opt<- par(mfrow=c(2,1))
RLE(Pset, main = "Relative Log Expression (RLE)", 
    names=rownames(pData(allAffyRaw)), las=2, cex.axis=0.6)
NUSE(Pset, main = "Normalized Unscaled Standard Errors (NUSE)", las=2, 
     names=rownames(pData(allAffyRaw)), las=2, cex.axis=0.6)

(#fig:plotPLM )Graficos de diagnóstico calculados a nivel de sondas PLM

par(opt)

4.3 Normalización de arrays de dos colores

La palabra describe las técnicas utilizadas para transformar adecuadamente los datos antes de que sean analizados. El objetivo es corregir diferencias sistemáticas entre muestras, en la misma o entre imágenes, lo que no representa una verdadera variación entre las muestras biológicas.

Estas diferencias sistemáticas pueden deberse, entre otras, a:

• Cambios en la tinción que producen sesgos la intensidad del .
• La ubicación en el array que puede afectar el proceso de lectura.
• Un problem en la placa de origen.
• La existencia de diferencias en la calidad de la impresión: pueden presentarse variaciones en los pins y el tiempo de impresión
• Camio en los parámetros de la digitalización (escaneo).

A veces puede ser difícil detectar estos problemas , aunque existen algunas formas de saber si es necesaria realizar una normalización. Aqui destacamos dos posibilidades:

1. Realizar una auto-hibridación. Si hibridamos una muestra con ella misma, las intensidades deberían ser las mismas en ambos canales. Cualquier desviación de esta igualdad, significa que hay un sesgo sistemático.
2. Detectar artefactos espaciales en la imagen o en la tinción de los gráficos de diagnóstico

4.3.1 Normalización global

Este método esta basado en un ajuste global, es decir en modificar todos los valores una cantidad , estimada de acuerdo a algún criterio. \[\begin{equation} \log_2 R/G \rightarrow \log_2 R/G-c=\log_2 R/(Gk) \end{equation}\] opciones para \(k\) o \(c= \log_2k\) son

\(c\)= mediana o media de log ratio para un conjunto concreto de genes o genes control o genes housekeeping.

La intensidad total de la normalización

\(k=\sum R_i/\sum G_i\)

4.3.2 Normalización dependiente de la intensidad

En este caso se realiza una modificación específica para cada valor. Esta modificación se obtiene como una función de la intensidad total del gen (\(c=c(A)\)). \[\begin{equation} \log_2 R/G \rightarrow \log_2 R/G-c(A)=\log_2 R/(Gk(A)) \end{equation}\]

Una posible estimación de esta función puede hacerse utilizando la función (LOcally WEighted Scatterplot Smoothing).

4.4 Resumen y normalización de microarrays de Affymetrix

En los arrays de Affymetrix, como en todos los tipos de microarrays, tras escanear la imagen se obtiene una serie de valores de intensidad de cada elemento del chip. En el caso de estos arrays sabemos que cada valor no corresponde a la expresión de un gen:

• Hay múltiples valores (sondas o “probes”) por cada gen, que originan un .
• Cada grupo de sondas consiste en múltiples pares de sondas, donde cada puede tener dos elementos: Un “perfect match” que coincide exactamente con el fragmento del gen al que corresponde la sonda Un “mismatch” que que coincide con el anterior salvo por el valor central que se ha sustituído por el nucleótido complementario. Estos "mismatches’ se introdujeron en los primeros arrays de affymetrix para tener una medida de hibridación no específica pero en las versiones más recientes se han abandonado.

El proceso que convierte las señales individuales en valores de expresión normalizados para cada gen consta de tres etapas:

1. Corrección del ruido de fondo o “background”
2. Normalización para hacer los valores comparables
3. “Sumarización”(Resumen) o concentración de los valores de cada grupo de sondas en un único valor de expresión absoluto normalizado para cada gen.

A menudo los tres pasos se denominan genérica -y erróneamente- “normalización”.

A diferencia de los chips de ADNc, aquí las medidas de expresión son absolutas (no se compara una condición contra otra) dado que cada chip se hibrida con un única muestra.

Hay muchos métodos para estimar la expresión (más de 30 publicados). Cada método contempla de forma explícita o implícita las tres formas de preprocesado: corrección del fondo, normalización y resumen.

Los principales métodos que consideraremos son:

• Microarray Suite (MAS). Método oficial de Affymetrix. Versiones 4.0 y 5.0
• dChip: Li and Wong. Método basado en modelos multichip.
• RMA (Bioconductor). Actualmente es el método “estándar”.

4.4.1 Métodos originales de Affymetrix

4.4.1.1 M.A.S. 4.0

Es el primer método introducida por Affymetrix. La corrección del fondo se realiza restando el “perfect match” del “mismatch” \[\begin{equation} E_j=PM_j-MM_j \end{equation}\] La normalización se realiza de forma global haciendo transformaciones de forma que la media de todo el chip sea la misma y la sumarización se basa en calcular el promedio de las diferencias absolutas ignorando los pares que se desvían más de \(3\sigma\) de \(\mu\). \[\begin{equation} Dif. Media=\frac{1}{|A|}\sum_{j \in A}(PM_j-MM_j) \end{equation}\] Los problemas que presenta estemétodo son:

• 1/3 de los MM son mayores que los PM
• Pueden aparecer valores MM negativos
• El uso de los MM añade ruido

4.4.1.2 M.A.S. 5.0

Los problemas que presentaba el método M.A.S. 4.0 llevaron a sustituirlo por otra variante, el M.A.S 5.0, llamado así por venir implementado en el software de affymetrix llamado “MicroArray Suite 5.0”.

Este método utiliza un estadístico robusto, , para corregir y ponderar el fondo y calcular (estimar) la señal. El biweight de Tukey \(T_{bi}\) pondera los valores por su distancia a la mediana, es decir, mide la tendencia central pero realiza un ajuste de outliers.

La lógica de este método reside en pensar que el valor de MM no siempre tiene sentido, (p.ej si MM \(>\) PM). Dado que esto sucede en ocasiones se realiza el cambio siguiente:

1. Se introduce el background específico de un conjunto de pruebas \(i\) de tamaño \(n\) basado en los pares de pruebas \(j\): \[\begin{equation} SB_i=T_{bi}(\log(PM_{i,j})-\log(MM_{i,j})): j= 1,\ldots,n \end{equation}\]
2. SB se utiliza para decidir como se ajusta el background

• Si es grande los datos suelen ser fiables
• Si es pequeño mejor basarse tan solo en PM

Este método no tan solo corrige el background sino que también permite normalizar y sumarizar. Para ello se introduce el “Mismatch Idealizado” (IM) que permite corregir la intensidad de las pruebas individuales. Este método también ha sido muy criticado:

• Se considera que no tiene mucho sentido promediar las pruebas entre arrays, pues éstos pueden tener características de hibridación intrínsecamente distintas.
• El método no mejora “aprendiendo” del funcionamiento entre arrays de las pruebas individuales.

4.4.2 El método RMA (Robust Multi-Array Average)

Para compensar algunas deficiencias de los primeros métodos de resumen y normalización de arrays de Affymetrix, Irizarry y sus colegas introdujeron en 2003 (~) un método basado en la modelización de las intensidades de las sondas que, en vez de basarse en las distintas sondas de un gen dentro de un mismo array se basa en los distintos valores de la misma sonda entre todos los arrays disponibles,

Esquemáticamente los pasos que realiza este método son:

1. Ajusta el ruído de fondo (background) basándose solo en los valores PM y utilizando un modelo estadístico complejo en el que combina la modelización de la señal mediante una distribución exponencial con la del ruído mediante una distribución normal.
2. Toma logaritmos base 2 de cada intensidad ajustada por el background.
3. Realiza una de los valores del paso 2 consistente en substituir cada valor individual por el que tendría la misma posición en la distribución empírica estimada sobre todas las muestras, es decir los promedios de las distribuciones de los valores ordenados de cada array (véase figura )
4. Estima las intensidades de cada gen separadamente para cada conjunto de sondas. Para ello realiza una técnica similar a una regresión robusta denominada “pulido de medianas” (median polish) sobre una matriz de datos que tiene los arrays en filas y los grupos de sondas en columnas.

Como resultado final de todos los pasos anteriores se obtiene la matriz con los datos sumarizados y normalizados. A pesar de no estar exento de críticas como la que afirma que este procedimiento “compacta” los valores reduciendo su variabilidad natural, este método se ha convertido en el estándar “de facto” actualmente por muchos usuarios de Bioconductor.

Figure 4.6: El método RMA incluye una normalización por cuantiles como la representada esquemáticamente en esta figura

4.5 Filtraje no específico

El filtraje no específico es recomendable para eliminar el ruido de fondo y limitar los ajustes posteriores a los necesarios. Los principales procesos de filtrado son:

• Eliminanción de los spots marcados como erróneos mediante flags y que son debidos a problemas en la hibridación o en el escaneo.
• Eliminanción de spots con señales muy bajas debido a problemas en el o a que no ha habido hibridación en ese spot (Filtraje pr ).
• Eliminación de genes que no presenten una variación significativa en su señal, entre distintas condiciones experimentales (Filtraje por variabilidad). Ante la duda,

El objetivo del filtraje es eliminar aquellos spots cuyas imágenes o señales sean erróneas por diferentes motivos, disminuyendo el ruido de fondo. Aunque existe controversia a su uso, prefiriendo el no filtrado a eliminar de forma no intencionado spots informativos.

5 Selección de genes diferencialmente expresados {# chapDEG}

5.1 Introducción

El motivo más habitual para el que se suelen utilizar microarrays es la búsqueda de genes cuya expresión cambia entre dos o más condiciones experimentales, por ejemplo a consecuencia de un tratamiento, una enfermedad u otras causas (distintos tiempos, distintas lineas celulares, …).

El problema consiste en identificar estos genes y suele denominarse o bien comparación de clases.

El problema de seleccionar genes diferencialmente expresados se traduce de manera casi inmediata al problema estadístico de comparar variables y, en años recientes, se han desarrollado un gran número de métodos estadísticos para resolverlo. La mayoría son extensiones de los métodos estadísticos clásicos, pruebas-t o análisis de la varianza, adaptados en uno u otro sentido para tener en cuenta las peculiaridades de los microarrays.

Aunque el problema de la selección de genes diferencialmente expresados puede relacionarse directamente con la realización de pruebas estadísticas, en el caso de los microarrays, el hecho de que haya dos tecnologías que miden la expresión de dos formas distintas hace que se deba diferenciar la metodología a emplear en cada caso. Los arrays de dos colores combinan dos muestras en un chip y generan una medida de expresión relativa. Esto hace que para comparar dos muestras de un mismo individuo sean la opción naturalmente más apropiada

En el caso de querer comparar muestras independientes de diferentes individuos los arrays de un color son la mejor opción. Evidentemente, lo más común será disponer de una sola técnica y tener que adaptar los análiaisis estadísticos a la misma.

Vamos a plantear un posible esquema de trabajo, para situar la mejor opción en cada caso:

• Situación 1: experimento con 5 individuos diabéticos y 5 no diabéticos, independientes entre si (muestras independientes)
  • Caso 1: Arrays de cDNA (2 colores): Utilizaríamos 5 arrays Diabético/Referencia y 5 arrays No diabético/Referencia
  • Caso 2: Arays de Affymetrix (1 color): Utilizaríamos 5 arrays de Diabético y 5 arrays de No diabético
• Situación 2: experimento con 6 individuos de los que se ha tomado una muestra de tejido sano y otra de tejido tumoral (muestras apareadas o dependientes)
  • Caso 3: Arrays de cDNA (2 colores): Utilizaríamos 6 arrays, uno por individuo, y en cada uno se realizaría la comparación Tejido Tumoral/Tejido sano.
  • Caso 4: Arays de Affymetrix (1 color): Utilizaríamos 12 arrays, 2 por individuo, 6 con muestras de Tejido Tumoral y 6 con muestras de Tejido sano.

5.2 Selección de genes diferencialmente expresados

5.2.1 Medidas naturales para comparar dos muestras

Recordemos que una vez se han hecho los experimentos con microarrays y obtenida la señal, los datos que se disponen son los logaritmos del valor detectado por el escaner.Esto hace que algunas operaciones que se realicen tengan en cuenta esta característica. Segun si la comparación a realizar se llevará a cabo con datos independientes (2 muestras, casos 1 y 2) o con datos dependientes (muestras apareadas, casos 3, 4) algunas medidas naturales o razonables para la comparación de expresiones son las siguientes:

• Para comparaciones directas, con expresiones relativas entre muestras apareadas o bien diferencias apareadas de expresiones absolutas:
  • log ratio promedio: \(\overline{R}=\frac 1n \sum_{i=1} R_i\)
  • t-test de una muestra \(\frac{\overline{R}}{SE}\), donde SE estima el error estándar del log ratio promedio
  • t-test robusto: Substituir en el anterior medidas robustas del error estándar
• Para comparaciones indirectas entre muestras independientes de expresiones relativas o absolutas:
  • Diferencia media \(\overline{R}_1-\overline{R}_2= \frac 1n_1 \sum_{i=1} ^{n_1}R_i- \frac 1n_2 \sum_{j=1} ^{n_2}R_j\)
  • t-test (clásico) de dos muestras \(\frac{\overline{R}_1-\overline{R}_2} {SE_{12}\sqrt{\frac 1n_1 +\frac 1n_2 }}\)
  • t-test robusto de dos muestras: Substituir en el anterior medidas robustas del error estándar

5.2.2 Un primer ejemplo

Consideremos la tabla siguiente que representa una matriz de expresión simplificada que contiene las expresiones relativas (por ejemplo entre tejido tumoral y sano del mismo individuo) de 5 genes en 6 muestras.

Podemos calcular las medidas descritas para el caso de una muestra para decidir si un gen está expresado o no lo está. Se discutirá más adelante como precisar esto pero de momento nos quedaremos con la idea de que si la medida escogida es (cercana a) cero el gen no está diferencialmente expresado y si es mayor o menor que cero si que lo está. Nos referimos a “cero” porque estamos hablando de logaritmos de razones: si la expresión es la misma en ambas condiciones el cociente es uno y su logaritmo es cero.

Vale la pena insistir en el concepto de expresión diferencial: no nos preocupa cual es la expresión del gen en una u otra muestra sinó si son distintas.

La tabla anterior sugiere que podría considerarse el gen A está diferencialmente expresado (promedio y t–test altos) mientras que el gen B o el D no lo están (promedio y test-t próximos a cero). Los genes C y D pueden llevar a conclusiones contradictoria según nos basemos en el promedio o el test t.

Si se observan los valores del test t del gen C se concluye que el gen no aparenta estar diferencialmente expresado. Si en cambio se observa su promedio parece que si que lo esté.

En el gen E pasa exctamente lo opuesto. La explicación de estas aparentes contradicciones se halla en el error estándar. En el gen C es muy elevado, debido a que el valor (20) es probablemente un “outlier”. En el gen D el error estándar es muy bajo por lo que, al encontrarse en el denominador del t-test aumenta artificialmente su valor.

5.2.3 Selección de genes usando tests estadísticos

Vamos a plantear la forma como se aborda este problema en el estudio de datos de microarrays. Dadas las características propias de este tipo de datos se consideran otras formas de estimar el error estándar que no sean tan sensibles a valores extremos o muy bajos.

Consideremos el gen \(g\). Si llamamos:

• \(R_g\) log-ratio medio observado.
• \(SE_g\) error estándar de \(R_g\) estimado a partir de los datos en el gen \(g\) .
• \(SE\) error estándar de \(R_g\) estimado a partir de los datos con la información de todos los genes.

Podemos considerar dos variantes para el test-t:

• Test-t global: se calcula en base a un único estimador de SE para todos los genes: \[t=R_g/SE,\]
• Test-t específico: Utiliza un estimador distinto del error estandar para cada gen: \[t=R_g/SE_g.\] Cada aproximación tiene sus pros y sus contras como muestra la tabla siguiente: \begin{center}

Test

Test-t Global

Test-t específico

\end{center}

En la práctica muchos métodos de selección de genes diferencialmente expresados han acabado buscando un compromiso entre ambas aproximaciones para lo que proponen o derivan fórmulas que de alguna forma ponderan o combinan dos estimaciones del error estándar: una basada en todos los genes y otra específica de cada gen. La tabla siguiente ilustra como algunos de los métodos más utilizados en la bibliografía incorporan esta idea.

Al hecho de incluir un coeficiente que tenga en cuenta la variabilidad de todos los genes en el array para estimar el error estándar de cada gen se le denomina moderación de la varianza (“variance shrinkage”) y es una de las aproximaciones en que existe cierto consenso (Allison et al. (2006)) acerca de que sirven para mejorar la selección de genes diferencialmente expresados.

5.2.3.1 Ejemplo de utilización del test-t

Como ejemplo utilizaremos el conjunto de datos celltypes y supondremos que disponemos ya de los datos normalizados y filtrados almacenados en un objeto expressionSet.

Este proceso es sencillo de hacer siguiendo los pasos del capiítulo anterior. Para simplificar el ejemplo una vez cargados y normalizados filtraremos los datos usando funcion nSFilter de la librería genefilter con sus opciones por defecto, lo que nos dejara un total de 10254 sondas en vez de las 45101 originales.

library(affy)
library(genefilter)
affyPath<- "datos/celltypes/celfiles"
cellTypesTargets<-  read.AnnotatedDataFrame("targets.txt", sep="\t",  path=affyPath)
cellTypesTargets$filename = file.path(affyPath, row.names(cellTypesTargets))
cellTypesRaw <- read.affybatch(cellTypesTargets$filename, phenoData=cellTypesTargets)
eset_rma <- rma(cellTypesRaw)
Filtered <- nsFilter(eset_rma)
eset_rma_filtered <- Filtered[["eset"]]
save(eset_rma, eset_rma_filtered, file="datos/celltypes/celltypes-normalized.rma.Rda")

Si consideramos que el campo treat del objeto contiene la información de los grupos a comparar (ratones estimulados con LPS frente a los que no han sido tratados), podemos realizar un primer análisis utilizando el test-t. Para ello utilizaremos el la función rowttests del paquete genefilter que realiza un test \(t\) sobre cada una de las filas (genes) de una matriz de expresión.

stopifnot(require(Biobase))
load (file="./datos/celltypes/celltypes-normalized.rma.Rda")
my.eset <- eset_rma_filtered
grupo_1 <- as.factor(pData(my.eset)$treat)
stopifnot(require(genefilter))
teststat <-rowttests(my.eset, "treat")
print(teststat[1:5,])

##                statistic         dm      p.value
## 1424378_at   -10.6739944 -1.1113660 8.715342e-07
## 1417675_a_at  20.7390631  0.9451576 1.504712e-09
## 1436530_at     9.2313983  1.4943922 3.291690e-06
## 1428603_at    -3.4882096 -0.4013432 5.840466e-03
## 1458888_at     0.7083721  0.2343463 4.948949e-01

Cuanto mayor sea el valor absoluto del estadístico \(t\) mayor es la probabilidad de que el gen esté diferencialmente expresado.

5.3 Significación estadística y expresión diferencial

Como hemos visto en la sección anterior dos medidas naturales para la selección de genes son el promedio de “log-ratios”, o la diferencia de promedios en el caso de muestras independientes, o el valor del estadístico de test (\(t\)-test) de una o dos muestras según si se trata de muestras apareadas o independientes respectivamente. Los primeros estudios de microarrays eran muy costosos y se hacían con pocas o incluso ninguna réplica por condición experimental. En estas situaciones la única forma fiable de detectar una diferencia de expresión era a traves del “log-ratio” o su diferencias.

Rápidamente se puso en evidencia que para poder obtener los genes que estaban realmente diferencialmente expresados era preciso disponer de un soporte estadístico que permitiera tener en cuenta la variación aleatoria existente entre muestras.

En la práctica esto se reduce a afirmar que si, además de la diferencia de expresión entre las condiciones experimentales, se lleva a cabo un test estadístico dispondremos de una medida objetiva, el p–valor que nos servirá para decidir qué genes se declaran diferencialmente expresados, a saber, aquellos en los que el p–valor del test sea inferior a un cierto umbral como 0.05 o 0.01.

Como es sabido, un test estadístico procede decidiendo rechazar la hipótesis nula si el p–valor es más pequeño que el nivel de significación del test.

Siguiendo con el ejemplo anterior podemos ordenar los resultados de los tests en base a los p–valores:

ranked <-teststat[order(teststat$p.value),]
print(ranked[1:5,])

##              statistic        dm      p.value
## 1449383_at   -48.61014 -2.044928 3.276616e-13
## 1451421_a_at -40.72173 -1.445182 1.909283e-12
## 1415929_at   -39.77036 -0.812879 2.415238e-12
## 1450826_a_at  37.62239  5.147992 4.193941e-12
## 1448303_at   -31.92365 -2.283765 2.140612e-11

Ahora podríamos seleccionar, por ejemplo los genes cuyo p–valor fuera inferior a 0.01

selectedTeststat <- ranked[ranked$p.value < 0.01,]

Esto deja un total de {r} dim(selectedTeststat)[1] con un p–valor inferior a 0.01

5.3.1 “Volcano plots”

Si se opta por computar los valores de significación (p–valores) de los genes, resulta interesante comparar el tamaño del cambio del nivel de significación estadístico. El “volcano plot” es una representación gráfica que permite ordenar los genes a lo largo de dos dimensiones, la biológica, representada por el “fold change” y la estadística representada por el logaritmo negativo del p–valor.

En la escala horizontal se representa el cambio entre los dos grupos (en escala logarítmica, de manera que la regulación positiva o negativa se representa de forma simétrica). En la escala vertical se representa el p–valor del test en una escala logarítmica negativa, de forma que los p–valores más pequeños aparecen mayores.

Así pues puede considerarse que el primer eje indica impacto biológico del cambio (a más efecto biológico mayor “fold-change”) y el segundo muestra la evidencia estadística, o la fiabilidad del cambio.

Como veremos más adelante en esta sección, el paquete limma permite realizar volcano-plots de forma sencilla a partir de los resultados de un anàlisis, pero para hacer un volcano plot tan sólo se necesita una lista de p-valores y los efectos (“fold-change”) asociados.

La figura (???)(fig:volcano0) muestra un “volcano-plot” para este ejemplo de los “celltypes” para la comparación “LPS” frente a “Medium”.

FC <- teststat$statistic
pVal<-teststat$p.value
Y <- -log(pVal)
plot(Y~FC)

También podemos adaptar alguna función “ad-hoc” como la siguiente, tomada de :

library(ggrepel)
# plot adding up all layers we have seen so far
volcanoP<- function (de,log2FoldChange,  pvalue, 
                     diffexpressed=NULL, col=NULL, delabel=NULL){
  ggplot(data=de, aes(x=log2FoldChange, y=-log10(pvalue), col=diffexpressed, label=delabel)) +
        geom_point() + 
        theme_minimal() +
        geom_text_repel() +
        scale_color_manual(values=c("blue", "black", "red")) +
        geom_vline(xintercept=c(-0.6, 0.6), col="red") +
        geom_hline(yintercept=-log10(0.05), col="red")
}
diffexpressed <- ifelse(abs(teststat$statistic)>10, TRUE, FALSE)
label <- rep(NA, length(diffexpressed))
label[diffexpressed] <- rownames(teststat)[diffexpressed]
volcanoP (de=teststat, log2FoldChange=teststat$statistic,  pvalue=teststat$p.value,
          diffexpressed = diffexpressed, delabel = label)

5.3.2 Potencia y tamaño muestral

Tal como se ha indicado más arriba, para realizar un suele controlarse la probabilidad de error de tipo I (de falsos positivos) y buscar, de entre los tests candidatos, aquellos que tengan una menor probabilidad de error de tipo II, o equivalentemente una mayor potencia. A partir de este planteamiento existe, en el contexto estadístico estándar, multitud de formas de determinar cual debe ser el tamaño muestral necesario para obtener una potencia dada fijados el tamaño de efecto (“fold-change”) y la probabilidad de error de tipo I.

En el caso de los microarrays se han desarrollado diversas fórmulas para realizar cálculos de este tipo, pero la compleja estructura de los datos microarrays hace que sean relativamente discutibles.

Simon (Simon et al. (2003)) sugiere la fórmula siguiente que es una generalización de las fórmulas clásicas para problemas de dos muestras:

El tamaño total requerido para detectar genes diferencialmente expresados en al menos una diferencia \(\delta\) con una probabilidad de error de tipo I (FP), \(\alpha\) y una probabilidad de error de tipo II (FN) \(1-\beta\) se calcula: \[ n=\frac{4(z_{\alpha/2}+z_\beta)^2}{(\delta/\sigma)^2}, \] donde \(z_\alpha\) y \(z_\beta\) son los percentiles 100\(\alpha/2\) y 100\(\beta\) de la distribución Normal \(N(0,1)\) y \(\sigma\) es la desviación estandar de un gen dentro de una clase (de un grupo). Obviamente \(\sigma\) es siempre desconocida por lo que, sin una muestra piloto con que estimarla el calculo e más imaginativo que realista.

Además de esto, el número de arrays usualmente recomendado queda lejos de la cantidad asequible para la mayor parte de los experimentos ((???),Tibshirani (2006)). Lo que muchos usuarios hacen a la práctica, es buscar un equilibrio entre los costes y la reproducibilidad y, de hecho, tienden a usar una cantidad fija de arrays tal como 3 o 5 sin consideraciones adicionales.

Por ejemplo si ponemos \(\alpha=0.001\), \(1-\beta=0.95\), \(\delta=1\) y estimamos \(\sigma\) entre todas las muestras el número de réplicas biológicas que necesitaremos será de 35.8.

5.4 El problema de la múltiplicidad de tests (“multiple testing”)

El análisis de microarrays se realiza en base gen a gen e involucra múltiples tests, miles probablemente. Esto significa que, a medida que crece el número de genes, la probabilidad de declarar erroneamente al menos un gen diferencialmente expresado va en aumento, y si no se realiza algún tipo de ajuste el número de falsos positivos será tanto más alto cuantos más genes estemos analizando.

Hay muchas formas de intentar controlar estas probabilidades de error y puede verse un excelente revisión en Dudoit Dudoit, Shaffer, and Boldrick (2003)).

De forma simplificada consideramos las dos aproximaciones más utilizadas.

Una posibilidad es mirar de controlar la probabilidad de obtener al menos un falso positivo o “Family-wise-error-rate (FWER)”. El más popular de estos métodos de control es la corrección de Bonferroni, consistente en multiplicar el p–valor por el número de tests realizados. La misma Dudoit y muchos otros autores han desarrollado variantes de los métodos clásicos de ajuste FWER usando por ejemplo tests de permutaciones.

El criterio FWER es quizás demasiado restrictivo dado que el control de los falsos positivos implica un considerable incremento de falsos negativos. En la práctica, sin embargo, muchos biólogos parecen estar dispuestos a aceptar que se produzcan algunos errores, siempre y cuando esto permita realizar descubrimientos. Por ejemplo un investigador debe considerar aceptable cierta pequeña proporción de errores (digamos del 10 al 20%) entre sus descubrimientos. En este caso, el investigador está expresando interés en controlar el porcentaje de falsos descubrimientos (FDR), es decir lo que es la proporción de falsos positivos sobre el total de genes inicialmente identificados como expresados diferencialmente. A diferencia del nivel de significación que queda determinado antes de examinar los datos, la FDR es una medida de confianza a posteriori ya que emplea información disponible en los datos para estimar las proporciones de falsos positivos que han tenido lugar. Si se obtiene una lista de los genes expresados diferencialmente en los que el FDR se controla hasta, digamos, el 20%, cabe esperar que el 20% de estos genes representen falsos positivos. Lo cual supone un enfoque menos restrictivo que controlar el FWER.

La decisión de controlar el FDR o el FWER depende de los objetivos del experimento. Si, por ejemplo, el objetivo es la es razonable permitir cierta cantidad de falsos positivos y es preferible seleccionar FDR. Si por el contrario se trabaja con una lista de un tamaño menor al deseado para verificar la expresión de ciertos genes específicos, entonces el FWER es el criterio apropiado.

El ejemplo siguiente muestra como se realiza el ajuste de p–valores usando el paquete multtest de Bioconductor para ajustar por los métodos de Bonferroni (FWER), Benjamini & Hochberg (FDR) o Benjamini & Yekutieli (BY, FDR).

stopifnot(require(multtest))
procs <- c("Bonferroni","BH", "BY")
adjPvalues <- mt.rawp2adjp(teststat$p.value, procs)
names(adjPvalues)

## [1] "adjp"    "index"   "h0.ABH"  "h0.TSBH"

ranked.adjusted<-cbind(ranked, adjPvalues$adjp)
head(ranked.adjusted)

##              statistic        dm      p.value         rawp   Bonferroni
## 1449383_at   -48.61014 -2.044928 3.276616e-13 3.276616e-13 3.359842e-09
## 1451421_a_at -40.72173 -1.445182 1.909283e-12 1.909283e-12 1.957778e-08
## 1415929_at   -39.77036 -0.812879 2.415238e-12 2.415238e-12 2.476585e-08
## 1450826_a_at  37.62239  5.147992 4.193941e-12 4.193941e-12 4.300467e-08
## 1448303_at   -31.92365 -2.283765 2.140612e-11 2.140612e-11 2.194983e-07
## 1418645_at   -28.74611 -4.126867 6.044555e-11 6.044555e-11 6.198087e-07
##                        BH           BY
## 1449383_at   3.359842e-09 3.296908e-08
## 1451421_a_at 8.255283e-09 8.100651e-08
## 1415929_at   8.255283e-09 8.100651e-08
## 1450826_a_at 1.075117e-08 1.054978e-07
## 1448303_at   4.389966e-08 4.307737e-07
## 1418645_at   1.033014e-07 1.013665e-06

Si seleccionamos los genes en base a su p–valor ajustado por ejemplo por le método de Banjamini y Yekutieli se obtienen los siguientes genes

selectedAdjusted<-ranked.adjusted[ranked.adjusted$BY<0.001,]
nrow(selectedAdjusted)

## [1] 420

5.5 Modelos lineales para la selección de genes: limma

En la sección anterior se ha discutido el uso del test \(t\) y sus extensiones para la selección de genes diferencialmente expresados en situaciones relativamente sencillas, es decir cuando sólo hay dos grupos.

En muchos estudios el número de grupos a considerar es más de dos y las fuentes de variabilidad pueden ser más de una, por ejemplo una puede ser el tratamiento pero otra puede ser la edad de los individuos o cualquier otra condición fijada por el investigador o derivadad de la heterogeneidad de las muestras.

En estas situaciones una aproximación razonable en problemas con una variable respuesta es el análisis de la varianza, discutido en el capítulo here. En esta sección se presenta una aproximación equivalente que de forma general se denomina el modelo lineal. Este método -que engloba el análisis d la varianza y la regresión- es uno de los más usados en estadística y ha sido popularizado en el campo de microarrays gracias a los trabajos de Gordon Smyth quien ha creado el paquete limma que se ha convertido en la herramiente más utilizada para el análisis de microarrays.

5.5.1 El modelo lineal general

El modelo lineal (ver por ejemplo Faraway, (???)) es un marco general para la modelización y el análisis de datos estadística.

EL método consiste en asumir una relación lineal entre los valores observados de una variable respuesta y las condiciones experimentales. A partir de aquí se obtienen estimadores para los parámetros del modelo y de sus errores estándar, y (con algunas condiciones extra) es posible hacer inferencia acerca del experimento.

La aplicación de modelos lineales puede ser visto como un proceso secuencial, con los siguientes pasos:

1. Especificar el diseño del experimento: qué muestras se asignan a qué condiciones.
2. (Re-)Escribir un modelo lineal para este diseño en forma de \(\mathbf{Y=X\beta+\varepsilon}\), donde \(\mathbf{X}\) se denomina la matriz de diseño.
3. Una vez que el modelo se ha especificado aplicar la teoría general de estimar los parámetros y los contrastes (comparaciones entre los valores de los parámetros).
4. Si se cumplen ciertas condiciones de validez para el modelo es posible realizar inferencia sobre los parámetros del modelo, es decir se pueden contrastar hipótesis sobre dichos parámetros.

El esquema anterior se puede aplicar a casi cualquier tipo de situación experimental. En la sección siguiente se presentan un par de ellas.

5.5.2 Ejemplos de situaciones modelizables linealmente

5.5.2.1 Ejemplo 1: Experimento “Swirl–Zebrafish”

Swirl es una mutación puntual que provoca defectos en la organización del embrión en desarrollo a lo largo de su eje dorsal-ventral. Como resultado, algunos tipos celulares se reducen y otros se expanden. Un objetivo de este trabajo fue identificar los genes con expresión alterada en el mutante Swirl en comparación con “wild zebrafish”.

1. El diseño experimental para este estudio fue el siguiente: \begin{center}

1
2
3
4

\end{center}

• Cada microarray contenía 8848 sondas de cDNA (genes o sequencias EST).
• Cuatro réplicas por array (slide): 2 juegos de pares de intercambio de color
• El cDNA del mutante swirl (\(S\)) se marca con Cy5 o Cy3 y el cDNA del “wild type” se marca con el otro

1. El modelo lineal derivado del diseño anterior fue:

• El parámetro de interés es: \(\alpha= \mathbf{E}\left (log\frac{S}{W}\right )\) .
• Las muestras 1 y 3 estan marcadas con : S (Verde=“Green”) y W (Rojo=“Red”), y las muestras 2 y 4 son “dye-swapped”.
• El modelo, \(\mathbf{y=X\alpha+\varepsilon}\), es: \[ \begin{array}{ccccc} y_1 &=&\alpha&+&\varepsilon _1 \\ y_2 &=&-\alpha&+&\varepsilon _2 \\ y_3 &=&\alpha&+&\varepsilon _3 \\ y_4 &=&-\alpha&+&\varepsilon _4\\ \end{array} <!-- \right.---> \Longrightarrow <!-- \left.---> \left( \begin{array}{r} y_1 \\ y_2 \\ y_3\\ y_4 \end{array} \right) = \underbrace{\left( \begin{array}{r} 1 \\ -1 \\ 1\\ -1 \end{array} \right)}_{\text{Matriz de Diseño}, \mathbf{X}} \alpha + \left( \begin{array}{r} \varepsilon_1 \\ \varepsilon_2 \\ \varepsilon_3\\ \varepsilon_4 \end{array} \right) \]

5.5.3 Ejemplo 2: Comparación de tres grupos

Los plásmidos IncHI codifican genes de resistencia múltiple a los antibióticos en S. enterica.

El plásmido R27 de la cepa salvaje es termosensible al transferirse.

Algunos fenotipos mutantes relacionados con hha y hns cromosómicos participan en diferentes procesos metabólicos de interés en la conjugación termoregulada.

El objetivo del experimento es encontrar qué genes se expresean diferencialmente en tres tipos de mutantes diferentes, \(M_1\), \(M_2\) y \(M_3\).

Este experimento debe ser planteado de forma diferente según el tipo de arrays (uno o dos colores) y qué comparaciones son las de mayor interés.

• Array de dos colores
  • A diseño de referencia: Hibridar cada Mutante (\(M_i\)) vs. Salvaje (“Wild type”) (\(W\)).
  • A diseño loop: Hibridar cada mutante el uno al otro en un doble “loop” que incluye (“dye-swapping”).
• Array de un color: hibridar mutantes y “wild types” separadamente.

" width=“50%”>

• Permite la comparación directa de Mutantes vs “Wild”.
• Número de parámetros a estimar es 3, relación intuitiva entre el número de parámetros y mutantes.
• Las comparaciones de Mutante vs Mutante son menos eficientes.

" width=“50%”>

• Permite la comparación directa de Mutantes vs Mutantes.
• El número de parámetros a estimar es 2. Menos intuitivo.
• Las comparaciones Mutante vs Mutante son más eficientes.

" width=“50%”>

• Permite la comparación directa de
  • Mutante vs “Wild”
  • Mutante vs Mutante
• El número de parámetros a estimar es 4.
• Todas las comparaciones se pueden hacer de forma eficiente.

Modelo, \(\mathbf{y=X\alpha+\varepsilon}\), y contrastes \(\mathbf{C'\beta}\)

• Parámetros del modelo: \[\alpha_1= \mathbf{E}\left (log\frac{M_1}{W}\right ),\ \alpha_2= \mathbf{E}\left (log\frac{M_2}{W}\right ),\ \alpha_3=\mathbf{E}\left (log\frac{M_3}{W}\right ).\]
• Contrastes: Comparaciones interesantes. \[\begin{array}{ccc} \beta_1 &=& \alpha_1-\alpha_2,\\ \beta_2 &=& \alpha_1-\alpha_3,\\ \beta_3 &=& \alpha_2-\alpha_3. \end{array}\]

\[ \left( \begin{array}{r} y_1 \\ y_2 \\ y_3 \\ y_4 \\ y_5 \\ y_6 \end{array} \right) = \underbrace{\left( \begin{array}{rrr} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \\ -1 & 0 & 0 \\ 0 & -1 & 0 \\ 0 & 0 & -1 \end{array} \right)}_{\text{Matriz de Diseño}, \mathbf{X}} \left( \begin{array}{c} \alpha_1 \\ \alpha_2 \\ \alpha_3 \end{array} \right) + \left( \begin{array}{r} \varepsilon_1 \\ \varepsilon_2 \\ \varepsilon_3\\ \varepsilon_4 \\ \varepsilon_5\\ \varepsilon_6 \end{array} \right) \] \[ \left( \begin{array}{r} \beta_1 \\ \beta_2 \\ \beta_3 \end{array} \right) = \underbrace{\left( \begin{array}{rrr} 1 & -1 & 0 \\ 1 & 0 & -1 \\ 0 & 1 & -1 \end{array} \right)}_{\text{Matriz de Contraste}, \mathbf{C}} \left( \begin{array}{c} \alpha_1 \\ \alpha_2 \\ \alpha_3 \end{array} \right). \]

Modelo, \(\mathbf{y=X\alpha+\varepsilon}\), y contrastes \(\mathbf{C'\beta}\)

• Parámetros del modelo: \[ \alpha_1= \mathbf{E}\left (log\frac{M_1}{M_2}\right ),\ \alpha_2= \mathbf{E}\left (log\frac{M_2}{M_3}\right ). \] \(\alpha_3\) no es necesaria: \(\log \left( \frac{M_1}{M_3}\right )=\log\left(\frac{M_1}{M_2}\right )-\log\left(\frac{M_2}{M_3}\right)\).

• Contrastes: Algunas de las comparaciones deseadas son precisamente los parámetros. \[\begin{array}{ccc} \beta_1 &=& \alpha_1,\\ \beta_2 &=& \alpha_2,\\ \beta_3 & =& \alpha_1+\alpha_2. \end{array}\]

\[ \left( \begin{array}{r} y_1 \\ y_2 \\ y_3 \\ y_4 \\ y_5 \\ y_6 \end{array} \right) = \underbrace{\left( \begin{array}{rr} 1 & 0 \\ 0 & 1 \\ 1 & -1 \\ -1 & 0 \\ 0 & -1 \\ -1 & 1 \end{array} \right)}_{\text{Matriz de Diseño}, \mathbf{X}} \left( \begin{array}{c} \alpha_1 \\ \alpha_2 \\ \end{array} \right) + \left( \begin{array}{r} \varepsilon_1 \\ \varepsilon_2 \\ \varepsilon_3\\ \varepsilon_4 \\ \varepsilon_5\\ \varepsilon_6 \end{array} \right) \]

\[ \left( \begin{array}{r} \beta_1 \\ \beta_2 \\ \beta_3 \end{array} \right) = \underbrace{\left( \begin{array}{rrr} 1 & 0 \\ 0 & 1 \\ 1 & +1 \end{array} \right)}_{\text{Matriz de Contraste}, \mathbf{C}} \left( \begin{array}{c} \alpha_1 \\ \alpha_2 \end{array} \right). \]

Modelo, \(\mathbf{y=X\alpha+\varepsilon}\), y contrastes \(\mathbf{C^{1'}\beta}\), \(\mathbf{C^{2'}\beta}\)

• Parámetros del modelo: \[\alpha_1= \mathbf{E} (log{M_1} ),\ \alpha_2= \mathbf{E} (log{M_2} ),\ \alpha_3=\mathbf{E} (log{M_3} ), \alpha_4=\mathbf{E} (log{W} ). \]
• Contrastes: Dos posibles conjuntos de comparaciones interesantes.

1. Comparación entre tipo de mutantes \((\mathbf{C^{1'}\beta})\) \[\begin{array}{ccc} \beta_1^1 &=& \alpha_1-\alpha_2,\\ \beta_2^1 &=& \alpha_3-\alpha_2,\\ \beta_3^1 &=& \alpha_2-\alpha_3. \end{array}\]

2. Comparación entre cada mutantes y el wild type \((\mathbf{C^{2'}\beta})\) \[\begin{array}{ccc} \beta_1^2 &=& \alpha_4-\alpha_1,\\ \beta_2^2 &=& \alpha_3-\alpha_1,\\ \beta_3^2 &=& \alpha_2-\alpha_1. \end{array}\]

\[ \left( \begin{array}{r} y_1 \\ y_2 \\ y_3 \\ y_4 \\ y_5 \\ y_6 \\ y_7 \\ y_8 \end{array} \right) = \underbrace{\left( \begin{array}{rrrr} 1 & 0 & 0 & 0\\ 1 & 0 & 0 & 0\\ 0 & 1 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \\ 0 & 0 & 0 & 1 \end{array} \right)}_{\text{Matriz de Diseño}, \mathbf{X}} \left( \begin{array}{c} \alpha_1 \\ \alpha_2 \\ \alpha_3 \\ \alpha_4 \end{array} \right) + \left( \begin{array}{r} \varepsilon_1 \\ \varepsilon_2 \\ \varepsilon_3\\ \varepsilon_4 \\ \varepsilon_5\\ \varepsilon_6\\ \varepsilon_7\\ \varepsilon_8 \end{array} \right) \]

\[ \left( \begin{array}{r} \beta_1^1 \\ \beta_2^1 \\ \beta_3^1 \end{array} \right) = \underbrace{\left( \begin{array}{rrrr} 1 & -1 & 0 & 0\\ 1 & 0 & -1 & 0\\ 0 & 1 & -1& 0 \end{array} \right)}_{\text{Matriz de Contraste}, \mathbf{C^1}} \left( \begin{array}{c} \alpha_1 \\ \alpha_2 \\ \alpha_3 \\ \alpha_4 \end{array} \right). \]

\[ \left( \begin{array}{r} \beta_1^2 \\ \beta_2^2 \\ \beta_3^2 \end{array} \right) = \underbrace{\left( \begin{array}{rrrr} 1 & 0 & 0 & -1\\ 0 & 1 & 0 & -1\\ 0 & 0 & 1 & -1 \end{array} \right)}_{\text{Matriz de Contraste}, \mathbf{C^2}} \left( \begin{array}{c} \alpha_1 \\ \alpha_2 \\ \alpha_3 \\ \alpha_4 \end{array} \right). \]

5.5.3.1 Ejemplo 3: Estudio de la influencia de las citoquinas en ratones viejos

El objetivo de este experimento es estudiar el efecto de las citoquinas sobre la estimulación de una sustancia (LPS) y ver como esta relación se ve afectada por la edad.

Se trata de un modelo de un un factor (tratamiento, asignable por el individuo) y un bloque (edad, no asignable, prefijada en cada ratón). En la práctica el modelo equivale al del análisis de la varianza de dos factores.

1. El diseño experimental para este estudio fue el siguiente: \begin{center}

% after : or …
1
2
3
4
5
6
7
8
9
10
11
12

\end{center}

• Se utilizáron microarrays de un color (Affymetrix).
• Cada condición se replicó tres veces.
• Las preguntas específicas a responder:
  • ?`Cual es el efecto del tratamiento en ratones viejos?
  • ?`Cual es el efecto del tratamiento en ratones jovenes?
  • ?`En que genes el efecto es diferente?.

1. En este caso se pueden considerar distintos modelo lineales derivados del hecho de que este experimento admite diferente parametrizaciones:

• Factores separados con 2 niveles cada uno para tratamiento (LPS/Med), Edad (Joven/Viejo) y su interacción: \[Y_{ijk}=\underbrace{\alpha_{i}}_{\text{Trat}}+\underbrace{\beta_{j}}_{\text{Edad}}+\underbrace{\gamma_{ij}}_{\text{Interacción}}+\varepsilon_{ijk},\, i=1,2,\, j=1,2,\, k=1,2,3\] Esta primera parametrización parece natural pero es más complicado confiar en ella para responder las preguntas propuestas.
• Un factor combinado con 4 niveles
  (LPS.Aged, Med.Aged, LPS.Young, Med.Young) \[ Y_{ij}=\alpha{i} +\varepsilon_{ij}, \quad i=1,...,4,\, j=1,2,3 . \] Esta parametrización parece más rígida pero se adapta mejor para responder a las preguntas planteadas.

Aquí se adopta la segunda parametrización.

• Parámetros del modelo: \[\begin{array}{ccc} \alpha_1&=& \mathbf{E} (log{LPS.Aged} ),\ \alpha_2= \mathbf{E} (log{Med.Aged} ),\\ \alpha_3&=&\mathbf{E} (log{LPS.Young} ),\, \alpha_4=\mathbf{E} (log{Med.Young} ). \end{array}\]

• Contrastes: Preguntas que interesa responder: \[\begin{array}{ccc} \beta_1^1 &=& \alpha_3-\alpha_1,\quad \mbox{Efecto del tratamiento en ratones viejos} \\ \beta_2^1 &=& \alpha_4-\alpha_2,\quad \mbox{Efecto del tratamiento en ratones jóvenes} \\ \beta_3^1 &=& (\alpha_3-\alpha_1)- (\alpha_2-\alpha_4),\quad \mbox{Interacción: diferencia entre efectos} \end{array}\]

• Modelo lineal:

Modelo, \(\mathbf{y=X\alpha+\varepsilon}\), y contrastes \(\mathbf{C'\beta}\)

$$ ( \[\begin{array}{r} y_1 \\ y_2 \\ y_3 \\ y_4 \\ y_5 \\ y_6 \\ y_7 \\ y_8 \\ y_9 \\ y_{10} \\ y_{11} \\ y_{12} \end{array}\] ) = _{, } ( \[\begin{array}{c} \alpha_1 \\ \alpha_2 \\ \alpha_3 \\ \alpha_4 \end{array}\] ) + ( \[\begin{array}{r} \varepsilon_1 \\ \varepsilon_2 \\ \varepsilon_3\\ \varepsilon_4 \\ \varepsilon_5\\ \varepsilon_6\\ \varepsilon_7\\ \varepsilon_8 \end{array}\]

) $$

\[ \left( \begin{array}{r} \beta_1^1 \\ \beta_2^1 \\ \beta_3^1 \end{array} \right) = \underbrace{\left( \begin{array}{rrrr} -1 & 0 & 1 & 0\\ 0 & -1 & 0 & 1\\ -1 & 1 & 1 & -1 \end{array} \right)}_{\text{Matriz de Contraste}, \mathbf{C}} \left( \begin{array}{c} \alpha_1 \\ \alpha_2 \\ \alpha_3 \\ \alpha_4 \end{array} \right) \]

5.5.4 Estimación e inferencia con el modelo lineal

Una vez se ha expresado el experimento como un modelo lineal: \[ \mathbf{E(y_g)=X\alpha_g},\quad \text{var}(y_g)=W_g\sigma_g, \] es posible usar la teoría estándar del modelo lineal (ver (???)) para obtener:

• Estimación de los parámetros: \(\hat \alpha_g( \approx \mathbf{\alpha})\).
• Desviación estándar de las estimaciones: \(\hat \sigma_g=s_g( \approx \sigma)\).
• Error estándar de las estimaciones:\(\widehat {\text{var}\hat \alpha_g}=V_g\,s_g^2\). Estas estimaciones son la base para realizar inferencia sobre \(\alpha\) i.e. test \(H_0:\ \alpha =0?\), basado en el hecho que: \[ t_{gj}=\frac{\alpha_{gj}}{s_g\sqrt{v_{gj}}}\sim \text{Distribución de Student}. \] De forma análoga pueden derivarse fórmulas para \(\alpha_1-\alpha_2\), es decir para decidir acerca de las comparaciones.

Los procedimientos de estimación y de inferencia no dependen de qué parametrización se ha adoptado, a pesar de que distintas parametrizaciones piedan dar lugar a distintos valores numéricos..

5.5.4.1 Fortaleza y debilidades del modelo lineal

El enfoque del modelo lineal es flexible y potente:

• Se puede adaptar a situaciones diferentes y complejas.
• Siempre produce buenas estimaciones (“BLUE”).
• Si las suposiciones son ciertas proporciona una base para la inferencia.

Por otro lado hay que tener en cuenta que, si las suposiciones no se cumplen, entonces las conclusiones deben tomarse con precaución.

Y lo que es peor, aún siendo válidas las suposiciones del modelo los resultados pueden verse afectados por el tamaño de la muestra de forma que, en muestras pequeñas donde pueden haber variaciones más grandes es fácil que se obtengan valores \(t\) no significativos o, por el contrario, excesivamente significativos (si la variación es muy reducida).

La metodología desarrollada por Smyth (Smyth, Michaud, and Scott (2005)) basada en los resultados de L"onsted & Speed (Speed (2003)) está dirigida a abordar cómo hacer frente a estas debilidades.

5.5.5 Modelos lineales para Microarrays

Smyth (Smyth, Michaud, and Scott (2005)) considera el problema de la identificación de genes que se expresan diferencialmente en las condiciones especificadas en el diseño de experimentos de microarrays. Como hemos dicho repetídamente la variabilidad de los valores de expresión difiere entre genes, pero la naturaleza paralela de la inferencia en microarrays sugiere la posibilidad de usar la información de todos los genes a la vez para mejorar la estimación de los parámetros, lo que puede llevar a una inferencia más fiable.

Básicamente lo que propone Smyth (Smyth, Michaud, and Scott (2005)) es una solución en tres pasos:

• Se plantea el problema como un model lineal con una componente bayesiana ya que se supone que los mismos parámetros a estimar son variables (no constantes) con distribuciones prior que se estimaran a partir de la información de todos los genes.
• A continuación se obtienen las estimaciones de los parámetros del modelo. La aproximación utilizada garantiza que estos estimadores tienen un comportamiento robusto incluso para pequeño número de arrays.
• Finalmente se calcula un “odd-ratio” que viene a ser la probabilidad de que un gen esté diferencialmente expresado frente a la de que no lo esté y se asocia este valor denominado estadístico \(B\) con un estadístico \(t\) moderado y su p–valor. \[ B=\log \frac{P[\text{Afectado}|M_{ij}]}{P[\text{No Afectado}|M_{ij}]}, \] gen=i \((i=1...N)\), réplica=j \((j=1,...,n)\).

El hecho de trabajar con logaritmos permite poner el punto de corte en el cero: A mayor valor positivo más probable es que el gen esté diferencialmente expresado. A mayor valor negativo, más probable es que no lo esté

5.5.6 Implementación y ejemplos

El paquete de Bioconductor limma al que hemos hecho referencia en los párrafos anteriores implementa el método de Smyth para la regularización de la varianza lo que lo ha convertido en muy popular entre los usuarios de microarrays

El código siguiente muestra como se crea la matriz del diseño y los contrastes para realizar el análisis mediante modelos lineales;

Empezamos por la matriz de diseño.

##                      Aged.LPS Aged.MED Young.LPS Young.MED
## Aged LPS 80L.CEL            1        0         0         0
## Aged LPS 86L.CEL            1        0         0         0
## Aged LPS 88L.CEL            1        0         0         0
## Aged Medium 81m.CEL         0        1         0         0
## Aged Medium 82m.CEL         0        1         0         0
## Aged Medium 84m.CEL         0        1         0         0
## Young LPS 75L.CEL           0        0         1         0
## Young LPS 76L.CEL           0        0         1         0
## Young LPS 77L.CEL           0        0         1         0
## Young Medium 71m.CEL        0        0         0         1
## Young Medium 72m.CEL        0        0         0         1
## Young Medium 73m.CEL        0        0         0         1
## attr(,"assign")
## [1] 1 1 1 1
## attr(,"contrasts")
## attr(,"contrasts")$lev
## [1] "contr.treatment"

Usando los nombres de las columnas de la matriz de diseño creamos los contrastes

require(limma)
cont.matrix <- makeContrasts (
      LPS.in.AGED=(Aged.LPS-Aged.MED),
      LPS.in.YOUNG=(Young.LPS-Young.MED),
      AGE=(Aged.MED-Young.MED),
      levels=design)
cont.matrix

##            Contrasts
## Levels      LPS.in.AGED LPS.in.YOUNG AGE
##   Aged.LPS            1            0   0
##   Aged.MED           -1            0   1
##   Young.LPS           0            1   0
##   Young.MED           0           -1  -1

Una vez definida la matriz de diseño y los contrastes podemos pasar a estimar el modelo, estimar los contrastes y realizar las pruebas de significación que nos indiquen, para cada gen y cada comparación, si puede considerarse diferencialmente expresado.

La penúltima instrucción ejecuta el proceso de regularización de la varianza utilizando modelos de Bayes empíricos para combinar la información de toda la matriz de datos y de cada gen individual y obtener estimaciones de error mejoradas.

El análisis proporciona los estadísticos de test habituales como Fold–change \(t\)-moderados o \(p\)-valores ajustados que se utilizan para ordenar los genes de mas a menos diferencialmente expresados.

A fin de controlar el porcentaje de falsos positivos que puedan resultar del alto numero de contrastes realizados simultaneamente los p–valores se ajustan de forma que tengamos control sobre la tasa de falsos positivos utilizando el metodo de Benjamini y Hochberg (Benjamini and Hochberg (1995)).

La funcion topTable genera para cada contraste una lista de genes ordenados de mas a menos diferencialmente expresados.

topTab_LPS.in.AGED <- topTable (fit.main, number=nrow(fit.main), coef="LPS.in.AGED", adjust="fdr")
topTab_LPS.in.YOUNG <- topTable (fit.main, number=nrow(fit.main), coef="LPS.in.YOUNG", adjust="fdr")
topTab_AGE  <- topTable (fit.main, number=nrow(fit.main) , coef="AGE", adjust="fdr")

Se puede visualizar los resultados mediante un volcano plot que, en este caso, representa en abscisas los cambios de expresión en escala logarítmica y en ordenadas el estadístico \(B\) en vez de el “menos logaritmo” del p-valor.

coefnum = 1
opt <- par(cex.lab = 0.7)
volcanoplot(fit.main, coef=coefnum, highlight=10, names=fit.main$ID,
            main=paste("Differentially expressed genes",colnames(cont.matrix)[coefnum], sep="\n"))
abline(v=c(-1,1))

par(opt)

6 Caso “Estrógeno”: Selección de genes diferencialmente expresados asociados con la resistencia a estrógenos.

6.1 Introducción

El análisis de datos de microarrays suele proceder, en general secuencialmente, a través de una serie de etapas tal y como se comentado en capítulos anteriores.

En cada una de las etapas pueden llevarse a cabo diversos procesos, con diversas variantes de métodos parecidos.A lo largo de la primera década del 2000 se ha desarrollado el proyecto Bioconductor que ha impulsado la creación de cientos de paquetes de R que implementan casi todas las capacidades posibles de análisis de microarrays. Aunque nadie puede conocer todos los paquetes con tan sólo conocer algunos de ellos es posible llevar a cabo potentes análisis con relativa facilidad.

El objetivo de este capítulo es ilustrar de forma breve pero completa como se hace un análisis de microarrays usando las facilidades que ofrecen Bioconductor y R.

Como se ha dicho, uno de los problemas en este tipo de estudios es que en cada etapa se puede proceder de varias formas lo que da lugar a un asfixiante número de posibilidades especialmente para el neófito. Con el fin de evitar este problema de momento se describe una sola de las opciones posibles para cada paso, lo que constituye un proceso “al estilo Bioconductor”

6.2 El estudio “Estrógeno”

El ejemplo analizado consiste en unos datos disponibles en forma de paquete (estrogen) en Bioconductor. Se trata de un estudio en el que se estudia la influencia del tratamiento con estrógeno y del tiempo transcurrido desde el tratamiento en la expresión de genes asociados con cancer de mama. En capítulos anteriores se han utilizado fragmentos de este ejemplo descrito en el capítulo here en el epígrafe corresponiente a los ejemplos (ejemplo here).

El diseño experimental para este estudio consiste en un diseño factorial de 2 factores “Estrogeno” (Presente o Ausente), “Tiempo” (10h o 48h). La documentación del paquete estrogen ofrece más detalles acerca del diseño y las variables utilizadas.

Si el paquete no se encuentra instalado puede instalarse con la instrucción BiocManager::install que se descarga de la web de Bioconductor.

if (!(require("estrogen", character.only=T))){
    BiocManager::install("estrogen")
    }

En lo que sigue supondremos que se ha llevado a cabo una instalación estándar de Bioconductor es decir se han ejecutado las instrucciones:

BiocManager::install()

6.2.1 Directorios y opciones de trabajo

Para facilitar supondremos que trabajamos en un directorio escogido por nosotros y cuya localización se asigna a la variable workingDir. Los datos se copiaran en un subdirectorio del anterior denominado “data” que se almacenará en la variable dataDir y los resultados se almacenarán en un directorio “results” cuyo nombre completo se almacenará en la variable resultsDir.

workingDir <-getwd()
if (!file.exists("datos")) system("mkdir datos")
if (!file.exists("datos/estrogen")) system("mkdir datos/estrogen")
if (!file.exists("results")) system("mkdir results")
dataDir <-file.path(workingDir, "datos/estrogen")
resultsDir <- file.path(workingDir, "results")
setwd(workingDir)

options(width=80)
options(digits=5)

6.3 Obtención y lectura de los datos

Para un análisis de datos de microarrays de Affymetrix se necesitan los archivos de imágenes escaneadas (“.CEL”) y un archivo en el que se asigne una condición experimental a cada archivo.

6.3.1 Los datos

Una ventaja de este ejemplo es que los datos se encuentran disponibles tras instalar el paquete estrogen por lo que se pueden copiar un directorio de trabajo o simplemente trabajar con ellos en el directorio original. El directorio en que se encuentran los archivos .CEL es:

library(estrogen)
estrogenDir <- system.file("extdata", package = "estrogen")
print(estrogenDir)

## [1] "/home/alex/R/x86_64-pc-linux-gnu-library/4.0/estrogen/extdata"

Para realizar el análisis es preciso copiar todos los archivos del directorio “extdata” a nuestro directorio de trabajo “data”.

ATENCION: Esta operación que depende de la instalación y del sistema operativo. Por ejemplo para copiarlos desde R en linux se escribirá:

system(paste ("cp ", estrogenDir,"/* ./data", sep=""))

Los archivos copiados son

• Los ocho archivos .CEL para el análisis
• Un archivo defectuoso para que se vea como debe salir un “array malo” en el control de calidad.
• Un archivo de covariables o “targets” que se usará al leer los archivos .CEL incorporar la información de las covariables en el objeto que se creará.

6.3.2 Lectura de los datos

El proceso de leer los datos puede parecer un poco extraño a primera vista pero la idea es simple:

En primer lugar creamos algunas estructuras de datos que contienen la información sobre las variables y - opcionalmente - sobre las anotaciones y el experimento y

A continuación nos basamos en estas estructuras para leer los datos y crear los objetos principales que se utilizarán para el análisis.

library(Biobase)
library(affy)
sampleInfo <- read.AnnotatedDataFrame(file.path(estrogenDir,"targLimma.txt"),
    header = TRUE, row.names = 1, sep="\t")
fileNames <- pData(sampleInfo)$FileName
rawData <- read.affybatch(filenames=file.path(estrogenDir,fileNames),
                          phenoData=sampleInfo)

En este ejemplo se incluye un archivo defectuoso con fines didácticos, llamado “badcel”. Para ilustrar las diferencias entre archivos correctos e incorrectos se puede leer en otro objeto.

library(affy)
sampleInfoWithBad <- read.AnnotatedDataFrame(file.path(dataDir,"phenoDataWithBad.txt"),
    header = TRUE, row.names = NULL, sep="\t")
fileNames <- pData(sampleInfoWithBad)$FileName
rawData.wrong <- read.affybatch(filenames=file.path(dataDir,fileNames),
                          phenoData=sampleInfoWithBad)

6.4 Exploración, Control de Calidad y Normalización

Tras leer los datos pasamos al preprocesado. Aunque puede interpretarse de distintas formas esta fase suele consistir en

1. Realizar algunos gráficos con los datos para hacerse una idea de como ha resultado el experimento.
2. Realizar un control de calidad más formal.
3. Normalizar y, en el caso de Affymetrix, resumir las expresiones

6.4.1 Exploración y visualización

La exploración previa puede hacerse paso a paso o en bloque si se utiliza algún paquete como arratQualityMetrics. En la ayuda de este paquete se encuentra una descripción de los análisis básicos que permite realizar.

El primer gráfico que suele hacerse és algun tipo de “densidad” que muestre la distribución de las señales en los datos “crudos”, es decir sin normalizar. Estos gráficos también permiten hacerse una idea de si las distribuciones de los distintos arrays son similares en forma y posición.

El gráfico de degradación –que no aparece en este caso, ya que daba problemeas al representarlo– permite hacerse una idea de como ha sido el proceso de hibridación de las muestras. Lineas paralelas sugieren una calidad similar.

##        low10A low10B  hi10A   hi10B   low48A   low48B   hi48A    hi48B
## slope  -0.103 -0.217 -0.129 -0.3810 -0.50400 -0.55300 -0.3510 -0.66700
## pvalue  0.550  0.194  0.394  0.0341  0.00482  0.00109  0.0464  0.00032

EL código para el gráfico de degradación sería:

El diagrama de cajas muestra, como el histograma, una idea de la distribución de los datos.

Finalmente un cluster jerárquico seguido de un dendrograma nos puede ayudar a hacernos una idea de si las muestras se agrupan por condiciones experimentales.

Si lo hacen es bueno, pero si no, no es necesariamente indicador de problemas, puesto que es un gráfico basado en todo los datos.

6.4.2 Control de calidad

Las exploraciones anteriores nos proporcionan una idea de como son los datos.

Se pueden realizar controles de calidad más estrictos como:

• Los controles de calidad estándar de Affymetrix, descritos en el paquete arrayQualityMetrics.

El paquete arrayQualityMetrics encapsula unos cuantos análisis, de forma que con una instrucción se pueden realizar todos los análisis y enviar la salida a un directorio.

if(!require(arrayQualityMetrics)) BiocManager::install("arrayQualityMetrics")
library(arrayQualityMetrics)
arrayQualityMetrics(rawData, outdir = "Informe_de_calidad_para_los_datos_del_caso_ESTROGENO")
arrayQualityMetrics(rawDataBad, outdir = "Informe_de_calidad_(2)_para_los_datos_del_caso_ESTROGENO")

Estos informes contienen gráficos, y explicaciones de como interpretarlos. Los gráficos llevan también indicadores para determinar si una muestra determinada puede tener algún tipo de problema.

En el código de ejemplo hemos realizado dos veces el análisis, una con los datos “buenos” y otra incluyendo el array defectuoso (“bad.cel”). Como puede verse en las tablas resumen, y en los directorios de resultados el array defectuoso es detectado por la mayoría de pruebas.

La figuras ?? y ?? presenta una tabla que resumen algunas de las pruebas (no todas) realizadas e indica si, para cada criterio, puede considerarse que un array dado es un outlier.

Para concluir con esta seccion merece la pena recordar que todos estos graficos son exploratorios. Raramente se descarta un array si solo un grafico lo sugiere.

6.4.3 Normalizacion y Filtraje

Una vez realizado el control de calidad se procede a normalizar los datos y sumarizarlos.

Hecho esto puede realizarse un filtraje no específico con el fin de eliminar genes que constituyen básicamente “ruído”, bien porque sus señales son muy bajas o bien porque apenas varían entre condiciones, por lo que no aportan nada a la selección de genes diferencialmente expresados.

La normalización puiede hacerse por distintos métodos (MAS5, VSN, RMA, GCRMA, …) En este ejemplo se utilizará el método RMA pero no se realizará filtraje alguno. Esto puede implicar quizás que para seleccionar genes diferencialmente expresados basándose en el ajuste de p-valores debamos utilizar criterios menos restrictivos que si hubiéramos filtrado, pero tiene la ventaja de eliminar un paso que, en el mejor de los casos, resulta controvertido.

El procesado mediante RMA implica un proceso en tres etapas:

• Corrección de fondo (el RMA hace precisamente esto).
• Normalización para hacer los valores de los arrays comparables.
• Summarización de las diversas sondas asociadas a cada grupo de sondas para dar un único valor.

library(affy)
eset_rma <- rma(rawData)

## Background correcting
## Normalizing
## Calculating Expression

Si se desea normalizar únicamente si los datos normalizados no estan disponibles puede utilizarse el código siguiente:

library(affy)
normalize <- T
if(normalize){
  eset_rma <- rma(rawData)
  save(eset_rma, file=file.path(dataDir,"estrogen-normalized.Rda"))
}else{
  load (file=file.path(dataDir,"estrogen-normalized.Rda"))
}

La normalización hace que los valores de los arrays sean comparables entre ellos, aunque los distintos métodos situan los valores en escalas distintas, por lo que lo que no resulta directamente comparable son los valores normalizados por distintos métodos.

Un boxplot de los valores normalizados muestra que los valores ya están en una escala en donde se pueden comparar.

Ahora bien, esto no debe considerarse sinónimo de buena calidad porque si se aplica la normalización por cuantiles el gráfico aparecerá así independientemente de otros factores que puedan afectar la calidad.

boxplot(eset_rma,main="Boxplot for RMA-normalized expression values ",
        names=sampleNames, cex.axis=0.7, col=info$grupo+1,las=2)

El código siguiente, que no se ejecuta muestra como se podrían aplicar distintos métodos para lego compararlos entre ellos.

eset_mas5 <- mas5(rawData)  # Uses expresso (MAS 5.0 method) much slower than RMA!
stopifnot(require(gcrma))
eset_gcrma <- gcrma(rawData) # The 'library(gcrma)' needs to be loaded first.
stopifnot(require(plier))
eset_plier <- justPlier(rawData, normalize=T) # The 'library(plier)' needs to be loaded first.
compara <-data.frame(RMA=exprs(eset_rma)[,1], MAS5 =exprs(eset_mas5)[,1],
                    GCRMA=exprs(eset_gcrma)[,1], PLIER =exprs(eset_plier)[,1])
pairs(compara)

6.4.3.1 Filtraje

El filtraje no específico permite eliminar los genes que varían poco entre condiciones o que deseamos quitar por otras razones como por ejemplo que no disponemos de anotación para ellos. La función nsFilter permite eliminar los genes que, o bien varían poco, o bien no se dispone de anotación para ellos.

Si al filtrar deseamos usar las anotaciones, o la falta de ellas, como criterio de filtraje debemos disponer del correspondiente paquete de anotaciones.

require(genefilter)
filtered <- nsFilter(eset_rma, require.entrez=TRUE,
         remove.dupEntrez=TRUE, var.func=IQR,
         var.cutoff=0.5, var.filter=TRUE,
         filterByQuantile=TRUE, feature.exclude="^AFFX")

La función nsFilter devuelve los valores filtrados en un objeto expressionSet y un informe de los resultados del filtraje.

names(filtered)

## [1] "eset"       "filter.log"

class(filtered$eset)

## [1] "ExpressionSet"
## attr(,"package")
## [1] "Biobase"

print(filtered$filter.log)

## $numDupsRemoved
## [1] 2856
## 
## $numLowVar
## [1] 4292
## 
## $numRemoved.ENTREZID
## [1] 1166
## 
## $feature.exclude
## [1] 19

eset_filtered <-filtered$eset

Podemos grabar el objeto eset_rma y los datos filtrados para su posterior uso.

save(eset_rma, eset_filtered, file=file.path(resultsDir, "estrogen-normalized.Rda"))

Despues del filtraje han quedado 4409 genes disponibles para analizar.

6.5 Selección de genes diferencialmente expresados

Como en las etapas anteriores la selección de genes diferencialmente expresados (GDE) puede basarse en distintas aproximaciones, desde la \(t\) de Student al programa SAM pasando por multitud de variantes.

En este ejemplo se aplicará la aproximación presentada por Smyth et al. (2004) basado en la utilización del combinada con un método para obtener una estimación mejorada de la varianza.

6.5.1 Análisis basado en modelos lineales

El capítulo 6 de tse manual o el manual del programa limma contienen explicaciones detalladas sobre construir un modelo lineal para este problema y como utilizarlo para seleccionar genes diferencialmente expresados.

6.5.1.1 Matriz de diseño

El primer paso para el análisis es crear la matriz de diseño.

La situación discutida en este ejemplo se puede modelizar de dos formas, tal como se discute en la presentación citada más arriba:

• Como un modelo de dos factores Estrogeno(Pres/Aus) Tiempo (10h/48h) con interacción.
• Como un modelo de un factor con cuatro niveles (Pres.10h/Pres.48h/Aus.10h/Aus.48h).

Tal como se describe en el manual de limma la segunda parametrización resulta a menudo más comoda, a pesar de parecer menos intuitiva, porque permite formular con más facilidad que la de dos factores las preguntas que típicamente interesan a los investigadores.

El modelo lineal para este estudio será:

\[ \left( \begin{array}{r} y_1 \\ y_2 \\ y_3 \\ y_4 \\ y_5 \\ y_6 \\ y_7 \\ y_8 \end{array} \right) = \underbrace{\left( \begin{array}{rrrr} 1 & 0 & 0 & 0\\ 1 & 0 & 0 & 0\\ 0 & 1 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \\ 0 & 0 & 0 & 1 \end{array} \right)}_{\mbox{Design Matrix}, \mathbf{X}} \left( \begin{array}{c} \alpha_1 \\ \alpha_2 \\ \alpha_3 \\ \alpha_4 \end{array} \right) + \left( \begin{array}{r} \varepsilon_1 \\ \varepsilon_2 \\ \varepsilon_3\\ \varepsilon_4 \\ \varepsilon_5\\ \varepsilon_6\\ \varepsilon_7\\ \varepsilon_8 \end{array} \right) \]

Los parámetros del modelos representan las cuatro combinaciones tiempo/estrogeno.

\[\begin{array}{ccc} \alpha_1&=& \mathbf{E} (log{Abs.10h} ),\\ \alpha_2&=& \mathbf{E} (log{Pres.10h} ),\\ \alpha_3&=& \mathbf{E} (log{Abs.48h} ),\\ \alpha_4&=& \mathbf{E} (log{Pres.48h} ). \end{array}\]

La matriz de diseño puede definirse manualmente o a partir de un factor creado específicamente para ello.

Manualmente, seria:

design.1<-matrix(
c(1,1,0,0,0,0,0,0,
  0,0,1,1,0,0,0,0,
  0,0,0,0,1,1,0,0,
  0,0,0,0,0,0,1,1),
nrow=8,
byrow=F)
colnames(design.1)<-c("neg10h", "est10h", "neg48h", "est48h")
rownames(design.1) <-  c("low10A", "low10B", "hi10A" , "hi10B",  "low48A", "low48B", "hi48A" , "hi48B")
print(design.1)

##        neg10h est10h neg48h est48h
## low10A      1      0      0      0
## low10B      1      0      0      0
## hi10A       0      1      0      0
## hi10B       0      1      0      0
## low48A      0      0      1      0
## low48B      0      0      1      0
## hi48A       0      0      0      1
## hi48B       0      0      0      1

Alternativamente puede crearse la matriz de diseño a partir de la información sobre las condiciones contenida en el phenoData, siempre que exista un campo adecuado para ello.

En este caso la columna Target se ha creado para utilizarla con esta finalidad. El objeto phenoData puede recrearse a partir del archivo original o extrayéndolo del objeto ExpresionSet que contiene los datos y las covariables.

##        neg10h est10h neg48h est48h
## low10A      1      0      0      0
## low10B      1      0      0      0
## hi10A       0      1      0      0
## hi10B       0      1      0      0
## low48A      0      0      1      0
## low48B      0      0      1      0
## hi48A       0      0      0      1
## hi48B       0      0      0      1
## attr(,"assign")
## [1] 1 1 1 1
## attr(,"contrasts")
## attr(,"contrasts")$lev
## [1] "contr.treatment"

Ambas matrices, design, design.1 resultan iguales.

6.5.1.2 Contrastes

Dado un modelo lineal definido a traves de una matriz de diseño pueden formularse las preguntas de interes como contrastes es decir comparaciones entre los parametros del modelo.

Cada parametrizacion distinta requerira de unos contrastes diferentes para las mismas preguntas, por lo que habitualmente se utilizara la parametrizacion que permita formular de forma mas clara las comparaciones de interes.

En este caso interesa estudiar

• Efecto del estrogeno al inicio del tratamiento
• Efecto del estrogeno al cabo de un tiempo del tratamiento
• Efecto del tiempo en ausencia de estrogeno

Esto se puede formular facilmente con la parametrizacion adoptada. \[\begin{array}{ccc} \beta_1^1 &=& \alpha_2-\alpha_1,\quad \mbox{Efecto del estrogeno pasadas 10 horas} \\ \beta_2^1 &=& \alpha_4-\alpha_3,\quad \mbox{Efecto del estrogeno pasadas 48 horas} \\ \beta_3^1 &=& \alpha_3-\alpha_1,\quad \mbox{Efecto del tiempo en ausencia de Estrogeno} \end{array}\]

cont.matrix <- makeContrasts (
      Estro10=(est10h-neg10h),
      Estro48=(est48h-neg48h),
      Tiempo=(neg48h-neg10h),
      levels=design)
cont.matrix

##         Contrasts
## Levels   Estro10 Estro48 Tiempo
##   neg10h      -1       0     -1
##   est10h       1       0      0
##   neg48h       0      -1      1
##   est48h       0       1      0

6.5.1.3 Estimación del modelo y selección de genes

Una vez definida la matriz de diseño y los contrastes podemos pasar a estimar el modelo, estimar los contrastes y realizar las pruebas de significación que nos indiquen, para cada gen y cada comparación, si puede considerarse diferencialmente expresado.

require(limma)
fit<-lmFit(eset_rma, design)
fit.main<-contrasts.fit(fit, cont.matrix)
fit.main<-eBayes(fit.main)

El método implementado en amplía el análisis tradicional utilizando modelos de Bayes empíricos para combinar la información de toda la matriz de datos y de cada gen individual y obtener estimaciones de error mejoradas.

El análisis proporciona los estadísticos de test habituales como Fold–change \(t\)-moderados o \(p\)-valores ajustados que se utilizan para ordenar los genes de mas a menos diferencialmente expresados.

A fin de controlar el porcentaje de falsos positivos que puedan resultar del alto numero de contrastes realizados simultaneamente los p–valores se ajustan de forma que tengamos control sobre la tasa de falsos positivos utilizando el metodo de Benjamini y Hochberg (Benjamini and Hochberg (1995)).

La funcion topTable genera para cada contraste una lista de genes ordenados de mas a menos diferencialmente expresados.

topTabEstro10 <- topTable (fit.main, number=nrow(fit.main), coef="Estro10", adjust="fdr")
topTabEstro48 <- topTable (fit.main, number=nrow(fit.main), coef="Estro48", adjust="fdr")
topTabTiempo  <- topTable (fit.main, number=nrow(fit.main) , coef="Tiempo", adjust="fdr")

Una forma de visualizar los resultados es mediante un volcano plot que representa en abscisas los cambios de expresión en escala logarítmica y en ordenadas el “menos logaritmo” del p-valor o alternativamente el estadístico \(B\).

6.5.2 Comparaciones múltiples

Cuando se realizan varias comparaciones a la vez puede resultar importante ver que genes cambian simultáneamente en más de una comparación. Si el número de comparaciones es alto también puede ser necesario realizar un ajuste de p-valores entre las comparaciones, distinto del realizado entre genes.

La función decidetests permite realizar ambas cosas. En este ejemplo no se ajustaran los p–valores entre comparaciones. Tan solo se seleccionaran los genes que cambian en una o más condiciones.

EL resultado del análisis es una tabla res que para cada gen y cada comparación contiene un 1 (si el gen esta sobre-expresado o “up” en esta condicion), un 0 (si no hay cambio significativo) o un -1 (si esta “down”-regulado).

Para resumir dicho análisis podemos contar qué filas tienen como mínimo una celda distinta de cero:

sum.res.rows<-apply(abs(res),1,sum)
res.selected<-res[sum.res.rows!=0,]
print(summary(res))

##        Estro10 Estro48 Tiempo
## Down        23      99     27
## NotSig   12512   12375  12591
## Up          90     151      7

Un diagrama de Venn permite visualizar la tabla anterior sin diferenciar entre genes “up” o “down” regulados.

\begin{figure}

vennDiagram (res.selected[,1:3], main="Genes in common #1", cex=0.9)

\end{figure}

7 Referencias

Alizadeh, A., M. B. Eisen, E. Davis, C. Ma, I. Lossos, A. Rosenwald, J. Boldrick, et al. 2000. “Distinct Types of Diffuse Large B–Cell Lymphoma Identified by Gene Expression Profiling.” Nature 403: 503–11.

Allison, David B, Xiangqin Cui, Grier P Page, and Mahyar Sabripour. 2006. “Microarray Data Analysis: From Disarray to Consolidation and Consensus.” Nat Rev Genet 7 (1): 55–65.

Benjamini, Y., and Y. Hochberg. 1995. “Controlling the False Discovery Rate: A Practical and Powerful Approach to Multiple Testing.” Journal of the Royal Statistical Society B 57: 289–300.

Chelvarajan, R. L., Y. Liu, D. Popa, M. L. Getchell, T. V. Getchell, A. J. Stromberg, and S. Bondada. 2006. “Molecular basis of age-associated cytokine dysregulation in LPS-stimulated macrophages.” Journal of Leukocyte Biology 79 (6): 1314.

Dudoit, S., J. P. Shaffer, and J. C. Boldrick. 2003. “Multiple Hypothesis Testing in Microarray Experiments.” Statistical Science 18: 71–103.

Golub, T. R., D. K. Slonim, P. Tamayo, C. Huard, M. Gaasenbeek, J. P. Mesirov, H. Coller, et al. 1999. “Molecular Classification of Cancer: Class Discovery and Class Prediction by Gene Expression Monitoring.” Science 286: 531–37.

Imbeaud, Sandrine, and Charles Auffray. 2005. “’The 39 Steps’ in Gene Expression Profiling: Critical Issues and Proposed Best Practices for Microarray Experiments.” Drug Discovery Today 10 (17): 1175–82. https://doi.org/10.1016/S1359-6446(05)03565-8.

Kerr, M K, and G A Churchill. 2001. “Experimental Design for Gene Expression Microarrays.” Biostatistics.

Moore, Shanna, Julia Vrebalov, Paxton Payton, and Jim Giovannoni. 2002. “Use of Genomics Tools to Isolate Key Ripening Genes and Analyse Fruit Maturation in Tomato.” Journal of Experimental Botany 53 (377): 2023–30. https://doi.org/10.1093/jxb/erf057.

Quackenbush, J. 2002. “Microarray Data Normalization and Transformation.” Nature Genet. 32: 496–501.

Scholtens, Denise, Alexander Miron, Faisal M. Merchant, Arden Miller, Penelope L. Miron, J. Dirk Iglehart, and Robert Gentleman. 2004. “Analyzing Factorial Designed Microarray Experiments.” J. Multivar. Anal. 90 (1): 19–43. https://doi.org/http://dx.doi.org/10.1016/j.jmva.2004.02.004.

Simon, Richard M., Edward L. Korn, Lisa M. McShane, Michael D. Radmacher, George W. Wright, and Yingdong Zhao. 2003. Design and Analysis of Dna Microarray Investigations. Springer-Verlag.

Smyth, Gordon K, Joëlle Michaud, and Hamish S Scott. 2005. “Use of Within-Array Replicate Spots for Assessing Differential Expression in Microarray Experiments.” Bioinformatics 21 (9): 2067–75.

Speed, T. 2003. Statistical Analysis of Gene Expression Data. Boca Raton, Fla.: Chapman & Hall/CRC.

Tibshirani, R. 2006. “A simple method for assessing sample sizes in microarray experiments.” BMC Bioinformatics 7 (1): 106.

van ’t Veer, Laura J, Hongyue Dai, Marc J van de Vijver, Yudong D He, Augustinus A M Hart, Mao Mao, Hans L Peterse, et al. 2002. “Gene Expression Profiling Predicts Clinical Outcome of Breast Cancer.” Nature 415 (6871): 530–6.

Yang, Y. H., M. J. Buckley, S. Dudoit, and T. P. Speed. 2002. “Comparison of Methods for Image Analysis on cDNA Microarray Data.” Journal of Computational and Graphical Statistic 11 (1).

Zhu, Qianqian, Jeffrey Miecznikowski, and Marc Halfon. 2010. “Preferred Analysis Methods for Affymetrix Genechips. II. An Expanded, Balanced, Wholly-Defined Spike-in Dataset.” BMC Bioinformatics 11 (1): 285. https://doi.org/10.1186/1471-2105-11-285.